################################################################################################################################################################################################################################################
# The Political Consequences of Ethnically Targeted Incarceration: Evidence from Japanese-American Internment During WWII
# Replication Code: Main Results by Camp and Partial SI
################################################################################################################################################################################################################################################

################################################################################################################################################################################################################################################
# Setup
################################################################################################################################################################################################################################################
# Packages
  require(haven)
  require(data.table)
  require(stargazer)
  require(ggplot2)
  require(nnet)
  require(RColorBrewer)
  require(maps)
  require(foreign)
  library(plyr)
  library(dplyr)
  library(magrittr)
  library(tidyr)
  library(stringr)
  library(xtable)
  library(ggjoy)
  library(RItools)
  library(lmtest)
  library(interplot)
  library(gridExtra)
  

# Functions
  splitit <- function(x,splitchar,n) sapply(strsplit(as.character(x), splitchar), "[", n)
  getmode <- function(v) {
    uniqv <- unique(v)
    uniqv[which.max(tabulate(match(v, uniqv)))]
  }
  as.numeric.factor <- function(x) {as.numeric(levels(x))[x]}
  
  cl  <- function(dat, fm, cluster){
    require(sandwich, quietly = TRUE)
    require(lmtest, quietly = TRUE)
    M <- length(unique(cluster))
    N <- length(cluster)
    K <- fm$rank
    dfc <- (M / (M - 1)) * ((N - 1) / (N - K))
    uj  <- apply(estfun(fm), 2, function(x) tapply(x, cluster, sum));
    vcovCL <- dfc * sandwich(fm, meat = crossprod(uj) / N)
    return(list(coeftest(fm, vcovCL), vcovCL))
    }

# Get p-value on F Statistic
  lmp <- function (modelobject) {
    if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
    f <- summary(modelobject)$fstatistic
    p <- pf(f[1],f[2],f[3],lower.tail=F)
    attributes(p) <- NULL
    return(p)
  }
  
################################################################################################################################################################################################################################################
# Directories
################################################################################################################################################################################################################################################
# This sets your working directory to the location of the source file. Note that this command is specific to RStudio. You'll need to change
# this if you're working directly in R or some other IDE (see https://stackoverflow.com/questions/13672720/r-command-for-setting-working-directory-to-source-file-location-in-rstudio)
  setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
  
################################################################################################################################################################################################################################################
# Load Data
################################################################################################################################################################################################################################################
# JARP Survey
  full <- readRDS("jarpfull_plus.RDS")

# WRA Data on Interned Population  
  wra <- read_dta("internment_raw.dta")
  wra <- data.frame(wra)

################################################################################################################################################################################################################################################
# Merge in Camp Information for Treatments
################################################################################################################################################################################################################################################
  camps <- sort(unique(full$camp[full$camp != "No Data" & full$camp != "Not evacuated or interned" & is.na(full$camp) == 0]))
  
  # Just looked counties up manually for these in order. Some places are in multiple counties, so the county listed is the one that contains most of the camp
  # Key to major camps: 
  # Gila River = Rivers, AZ in Pinal County, AZ
  # Granada = Amache, CO in Prowers County, CO
  # Heart Mountain = Heart Mountain, Wyoming in Park County
  # Jerome = Denson, AR mostly in Drew but partly in Chicot county
  # Manzanar = Manzanar, CA in Inyo County
  # Minidoka = Hunt, ID in Jerome County, Idaho
  # Poston = Poston, AZ in Yuma County
  # Rohwer = McGehee, AR in Desha County Arkansas
  # Topaz = Topaz, UT in Utah
  # Newell, CA is Tule Lake, which is mostly Modoc but partly in Siskiyou county
  # Camp Livingston, LA is on the Rapides Parish and Grant Parish Lines
  
  # Merge in County Names for Various Camps
  camps <- data.frame(camps = camps, V1 = cbind(splitit(camps, ",", 1)))
  camps$counties <- c("Prowers", "Burleigh", "Zavala", "Drew", "Comanche", "Park", "Jerome", "Rapides", "Hidalgo", "Inyo", "Desha", "Missoula", "Modoc", "Yuma", "Pinal", "Santa Fe", "Millard")
  camps$stateabb <- splitit(camps$camps, ", ", 2)
  camps <- merge(camps, data.frame(cbind(state.abb, state.name)), by.x = "stateabb", by.y = "state.abb")
  camps <- camps[, c("camps", "counties", "state.name")]
  colnames(camps) <- c("camp", "county", "state")
  camps$county <- toupper(camps$county) ; camps$state <- toupper(camps$state)
  
  # Identify Additional Camp Treatments
  # Watch Towers, Military Police Buildings, Strikes (see Confinement and Ethnicity book in References)
  # towers = watch towers in camp; bldgs = military buildings on premises; strikes = strikes or work stoppages;
  # force = camp personnel use force or violence against internees; violence = violence among internees or between internees and civilian population;
  # peakpop = peak population (Table 3.2 in Confinement and Ethnicity)
  # Other sources: http://newmexicohistory.org/places/lordsburg-internment-pow-camp (Lordsburg), http://encyclopedia.densho.org/Fort_Sill_(detention_facility)/
  # 
  severity <- data.frame(camp = c("Rivers, AZ", "Amache, CO", "Heart Mountain, WY", "Denson, AR", "Manzanar, CA", "Hunt, ID",
                                  "Poston, AZ", "McGehee, AR", "Topaz, UT", "Newell, CA", "Fort Sill, OK", "Lordsburg, NM", "Bismarck, ND",
                                  "Crystal City, TX", "Livingston, LA"),
                         campname = c("Gila River", "Amache", "Heart Mountain", "Jerome", "Manzanar", "Minidoka",
                                      "Poston", "Rohwer", "Topaz", "Tule Lake", "Fort Sill", "Lordsburg", "Fort Lincoln",
                                      "Crystal City", "Livingston"),
                         towers = c(1, 6, 9, 7, 8, 8, 0, 8, 7, 19, NA, NA, NA, NA, NA),
                         bldgs = c(15, 15, 19, 12, 13, 15, 12, 12, NA, 63, NA, NA, NA, NA, NA),
                         strikes = c(0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1),
                         force = c(0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0),
                         violence = c(0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),
                         peakpop = c(13348, 7318, 10767, 8497, 10046, 9397, 17814, 8475, 8130, 18789, 707, 2500, 1518, 4000, 1123))
  
  severity$militarism <- (severity$towers) / (severity$peakpop / 1000)
  
  severity$stateview <- ifelse(severity$strikes == 1 | severity$force == 1, 1, 0)
  
  severity$any <- ifelse(severity$strikes == 1 | severity$force == 1 | severity$strikes == 1, 1, 0)
  
  camps <- merge(camps, severity, by = "camp")
  
  full <- merge(full, camps, by = "camp", all.x = T)
  
  full$interned.numeric <- ifelse(full$interned %in% c("Yes", "Relocated (evacuated)", "1-2 years", "1 year or less", "2-3 years",
                                                       "3 years or more", "If both 2 and 3", "Interned"), "Yes",
                                  ifelse(full$interned %in% c("None", "No", "Not evacuated or interned", "Stated code 3 (detention; claimed relocation but only detained", "Voluntary evacuation"), "No", NA))  
  
  full$interned.length.numeric <- ifelse(full$interned.length %in% c("Less than one year", "1 year or less", "2 to 6 months", "Under 2 months", "6 to 12 months"), 6,
                                         ifelse(full$interned.length %in% c("1 but less than 2 years", "1-2 years", "1 to 2 years"), 18,
                                                ifelse(full$interned.length %in% c("2 but less than 3 years", "2-3 years", "2 to 3 years"), 30,
                                                       ifelse(full$interned.length %in% c("3 but less than 4 years", "3 to 4 years", "2 to 3 years", "3 years or more"), 42,
                                                              ifelse(full$interned.length %in% c("4 but less than 5 years", "4 years or more"), 54,
                                                                     ifelse(full$interned.length %in% c("5 years or more"), 66, NA))))))
  
  # Summary Statistics for Internment Length
    length.summary <- data.frame(table(full$interned.length.numeric))
    length.summary$prop <- round(length.summary$Freq / sum(length.summary$Freq) * 100, 1)
    length.summary$Var1 <- c("Less than 1 Year", "1 - 2 Years", "2 - 3 Years", "3 - 4 Years", "4 - 5 Years", "5 + Years")
    colnames(length.summary) <- c("Internment Length", "N", "Percentage of Interned Sample")
    
  # Manuscript Table 3  
    stargazer(length.summary, summary = F, rownames = F, digits = 2, title = "Internment Length - Summary Statistics",
              out = "ms_table3.tex")
    
################################################################################################################################################################################################################################################
# Recode Outcomes & Covariates
################################################################################################################################################################################################################################################
# Demographics
  # Age
    full$age <- as.numeric(full$age)
  
  # One respondent has age of 0, which has to be an error in the data, so set to NA
    full$age[full$respid == "0845" & full$birth.order == "60"] <- NA

# Generate indicator for pre-internment location across generations
  full$i.precamploc <- as.character(full$i.precamploc) ; full$n.precamp <- as.character(full$n.precamp) ; full$s.precamp <- as.character(full$s.precamp)
  
  full$precamp <- ifelse(full$generation == "Issei", full$i.precamploc,
                         ifelse(full$generation == "Nisei", full$n.precamp,
                                ifelse(full$generation == "Sansei", full$s.precamp, NA)))
  
  # Generate indicator for farm vs. nonfarm occupation. For Issei, this is going to be the head of household's occupation from 1932-1941;
  # for Nisei it will be the respondent's own primary occupation from 1932-1941
  # Convert to numeric  
    full$i.pthhjob <- as.numeric(full$i.pthhjob) ; full$n.ptjob <- as.numeric(full$n.ptjob)
  
  # Generate single indicator because these are both coded via note 3
    full$ptjob <- ifelse(full$generation == "Issei", full$i.pthhjob,
                         ifelse(full$generation == "Nisei", full$n.ptjob, NA))
    
  # Flag farm occupations
    full$ptfarm <- ifelse(full$ptjob %in% c(222, 901, 905, 981, 982, 983, 984, 985, 986, 987, 988, 989), 1, 0)
    
# Define Treatment Status    
  full$anyfam <- ifelse(full$internedinfam > 0, 1, 0)
  
  full$treatmentstatus <- ifelse(full$anyfam == 1 & full$interned.numeric == "Yes", "self.fam",
                                 ifelse(full$anyfam == 0 & full$interned.numeric == "Yes", "self.only",
                                        ifelse(full$anyfam == 1 & full$interned.numeric == "No", "fam.only",
                                               ifelse(full$anyfam == 0 & full$interned.numeric == "No", "control", NA))))
  
# Create Dummies for California and being a minor
  full$cadummy <- ifelse(substr(full$precamp, 1, 2) %in% c("93", "94", "95"), 1, 0)
  full$minordummy <- ifelse(full$age - 25 < 18, 1, 0)
  
# Create Single Birthplace Variable for Sansei and Nisei  
  full$birthplace <- ifelse(full$generation == "Issei", NA,
                            ifelse(full$generation == "Nisei", full$n.bplace,
                                   ifelse(full$generation == "Sansei", full$s.bplace, NA)))

################################################################################################################################################################################################################################################
# Results: Level of Interest in U.S. Politics
################################################################################################################################################################################################################################################
# Recode Political Interest. This is asked differently of the three generations; the Issei just say Y/N and the S and N give a range of responses
  interest <- unique(c(full$i.interestus, full$n.interestus, full$s.interestus))
  interestcodes <- rep(NA, length(interest))

  interest <- data.frame(cbind(interest, interestcodes))
  interest$interestcodes <- ifelse(interest$interest == "No" | interest$interest == "No interest at all", 0,
                                   ifelse(interest$interest == "Only a little" | interest$interest == "Yes", 1,
                                          ifelse(interest$interest == "A fair amount", 2,
                                                 ifelse(interest$interest == "A great deal", 3, NA))))
  
  interest$interestcodes.binary <- ifelse(interest$interest %in% c("No", "No interest at all"), 0,
                                          ifelse(is.na(interest$interest), NA, 1))
  
  interest$interest <- as.character(interest$interest)

  full$ptxinterest <- ifelse(full$generation == "Issei", full$i.interestus,
                             ifelse(full$generation == "Nisei", full$n.interestus, 
                                    ifelse(full$generation == "Sansei", full$s.interestus, NA)))

  full <- merge(full, interest, by.x = "ptxinterest", by.y = "interest")
  
  full$exposure <- ifelse(full$treatmentstatus %in% c("self.only", "self.fam"), "self.interned", 
                          ifelse(full$treatmentstatus %in% c("fam.only"), "fam.only", "control"))
  
# Run Regressions
  # Violence
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
  
    # Model  
      interest.violence <- lm(interestcodes ~ violence*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.violence.ses  <- cl(comp, interest.violence, comp$camp)
      
    # Violence-Youth
      # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("interestcodes", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
      # Model  
      interest.violence.youth <- lm(interestcodes ~ violence*exposure + precamp + gender + age + generation, data = comp)  
      
      # Clustered Standard Errors 
      interest.violence.youth.ses  <- cl(comp, interest.violence.youth, comp$camp)
      
  # Violence - Age Interaction
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("interestcodes", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.violence.ageinter <- lm(interestcodes ~ violence * age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.violence.ageinter.ses  <- cl(comp, interest.violence.ageinter, comp$camp)    
      
  # Violence - Minors
      # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("interestcodes", "violence", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
      comp <- comp[complete.cases(comp),]
      
      # Model  
      interest.violence.minors <- lm(interestcodes ~ violence * minordummy + age + precamp + gender + generation, data = comp)  
      
      # Clustered Standard Errors 
      interest.violence.minors.ses  <- cl(comp, interest.violence.minors, comp$camp)        
  
  # Violence - Birthplace Controls
      # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "violence", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
      # Model  
      interest.violence.birthplace <- lm(interestcodes ~ violence*exposure + birthplace + gender + age + generation, data = comp)  
      
      # Clustered Standard Errors 
      interest.violence.birthplace.ses  <- cl(comp, interest.violence.birthplace, comp$camp)
      
  # Violence - LA Only
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "violence", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      comp <- comp[comp$birthplace == "933",]
      
    # Model  
      interest.violence.la <- lm(interestcodes ~ violence*exposure + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.violence.la.ses  <- cl(comp, interest.violence.la, comp$camp)
      
  # Force
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.force <- lm(interestcodes ~ force*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.force.ses  <- cl(comp, interest.force, comp$camp)
      
    # Force - Youth
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("interestcodes", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        interest.force.youth <- lm(interestcodes ~ force*exposure + precamp + gender + age + generation, data = comp)  
      
      # Clustered Standard Errors 
        interest.force.youth.ses  <- cl(comp, interest.force.youth, comp$camp)
    
    
    # Force - Age Interaction    
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("self.interned"), c("interestcodes", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        interest.force.ageinter <- lm(interestcodes ~ force * age + precamp + gender + generation, data = comp)  
        
      # Clustered Standard Errors 
        interest.force.ageinter.ses  <- cl(comp, interest.force.ageinter, comp$camp)    
        
    # Force - Minors
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("self.interned"), c("interestcodes", "force", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        interest.force.minors <- lm(interestcodes ~ force * minordummy + age + precamp + gender + generation, data = comp)  
        
      # Clustered Standard Errors 
        interest.force.minors.ses  <- cl(comp, interest.force.minors, comp$camp)       
        
    # Force - Birthplace Controls
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "force", "birthplace", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        interest.force.birthplace <- lm(interestcodes ~ force*exposure + birthplace + gender + age + generation, data = comp)  
        
      # Clustered Standard Errors 
        interest.force.birthplace.ses  <- cl(comp, interest.force.birthplace, comp$camp)
        
    # Force - LA Only
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "force", "birthplace", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        comp <- comp[comp$birthplace == "933",]
        
      # Model  
        interest.force.la <- lm(interestcodes ~ force*exposure + gender + age + generation, data = comp)  
        
      # Clustered Standard Errors 
        interest.force.la.ses  <- cl(comp, interest.force.la, comp$camp)    
        
  # Strikes 
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.strikes <- lm(interestcodes ~ strikes*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.strikes.ses  <- cl(comp, interest.strikes, comp$camp)
      
  # Strikes - Youth
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only")  & full$age - 22 < 17, c("interestcodes", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.strikes.youth <- lm(interestcodes ~ strikes*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.strikes.youth.ses  <- cl(comp, interest.strikes.youth, comp$camp)
      
  # Strikes - Age Interaction
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("interestcodes", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.strikes.ageinter <- lm(interestcodes ~ strikes*age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.strikes.ageinter.ses  <- cl(comp, interest.strikes.ageinter, comp$camp)    
      
  # Strikes - Minors
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("interestcodes", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.strikes.minors <- lm(interestcodes ~ strikes*minordummy + age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.strikes.minors.ses  <- cl(comp, interest.strikes.minors, comp$camp)   
      
  # Strikes - Birthplace Controls
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "strikes", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.strikes.birthplace <- lm(interestcodes ~ strikes*exposure + birthplace + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.strikes.birthplace.ses  <- cl(comp, interest.strikes.birthplace, comp$camp)
      
  # Strikes - LA Only
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "strikes", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      comp <- comp[comp$birthplace == "933",]
      
    # Model  
      interest.strikes.la <- lm(interestcodes ~ strikes*exposure + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.strikes.la.ses  <- cl(comp, interest.strikes.la, comp$camp)
      
      
  # Militarism
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.militarism <- lm(interestcodes ~ militarism*exposure + precamp + gender + age + generation, data = comp)  
    
    # Clustered Standard Errors 
      interest.militarism.ses  <- cl(comp, interest.militarism, comp$camp)
      
  # Militarism - Youth
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("interestcodes", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.militarism.youth <- lm(interestcodes ~ militarism*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.militarism.youth.ses  <- cl(comp, interest.militarism.youth, comp$camp)
      
  # Militarism - Age Interaction
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("interestcodes", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.militarism.ageinter <- lm(interestcodes ~ militarism*age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.militarism.ageinter.ses  <- cl(comp, interest.militarism.ageinter, comp$camp)    
      
  # Militarism - Minors
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("interestcodes", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.militarism.minors <- lm(interestcodes ~ militarism*minordummy + age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.militarism.minors.ses  <- cl(comp, interest.militarism.minors, comp$camp)      
      
  # Militarism - Birthplace Controls
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "militarism", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      interest.militarism.birthplace <- lm(interestcodes ~ militarism*exposure + birthplace + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.militarism.birthplace.ses  <- cl(comp, interest.militarism.birthplace, comp$camp)
      
  # Militarism - LA Only
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("interestcodes", "militarism", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      comp <- comp[comp$birthplace == "933",]
      
    # Model  
      interest.militarism.la <- lm(interestcodes ~ militarism*exposure + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      interest.militarism.la.ses  <- cl(comp, interest.militarism.la, comp$camp)    
      
################################################################################################################################################################################################################################################
# Results: Faith in Government
################################################################################################################################################################################################################################################
# This question is phrased "Most people in government are not really interested in the problems of the average man". So agreeing here means you don't think much of government.
  full$govhelp <- ifelse(full$generation == "Nisei", full$n.govhelp,
                         ifelse(full$generation == "Sansei", full$s.govhelp, NA))

# Numeric coding
  full$govhelp <- ifelse(full$govhelp == "Agree", 1,
                         ifelse(full$govhelp == "Disagree", 0, NA))

# Run Regressions
  # Violence  
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
    
    # Model  
      faith.violence <- lm(govhelp ~ violence*exposure + precamp + age + gender + generation, data = comp)
    
    # Clustered Standard Errors 
      faith.violence.ses  <- cl(comp, faith.violence, comp$camp)
      
  # Violence - Youth
      comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("govhelp", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.violence.youth <- lm(govhelp ~ violence*exposure + precamp + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.violence.youth.ses  <- cl(comp, faith.violence.youth, comp$camp)
  
  # Violence - Age Interaction 
      comp <- full[full$exposure %in% c("self.interned"), c("govhelp", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.violence.ageinter <- lm(govhelp ~ violence*age + precamp + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.violence.ageinter.ses  <- cl(comp, faith.violence.ageinter, comp$camp) 
      
  # Violence - Minors
      comp <- full[full$exposure %in% c("self.interned"), c("govhelp", "violence", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.violence.minors <- lm(govhelp ~ violence*minordummy + age + precamp + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.violence.minors.ses  <- cl(comp, faith.violence.minors, comp$camp)    
      
  # Violence - Birthplace Controls  
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "violence", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.violence.birthplace <- lm(govhelp ~ violence*exposure + birthplace + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.violence.birthplace.ses  <- cl(comp, faith.violence.birthplace, comp$camp)  
      
  # Violence - LA Only  
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "violence", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      comp <- comp[comp$birthplace == "933",]
      
    # Model  
      faith.violence.la <- lm(govhelp ~ violence*exposure + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.violence.la.ses  <- cl(comp, faith.violence.la, comp$camp)      
      
  # Force
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.force <- lm(govhelp ~ force*exposure + precamp + age + gender + generation, data = comp)
    
    # Clustered Standard Errors 
      faith.force.ses  <- cl(comp, faith.force, comp$camp)
      
  # Force - Youth
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("govhelp", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.force.youth <- lm(govhelp ~ force*exposure + precamp + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.force.youth.ses  <- cl(comp, faith.force.youth, comp$camp)
      
  # Force - Age Interaction
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("govhelp", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.force.ageinter <- lm(govhelp ~ force*age + precamp + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.force.ageinter.ses  <- cl(comp, faith.force.ageinter, comp$camp)    
      
  # Force - Minors
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("govhelp", "force", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.force.minors <- lm(govhelp ~ force*minordummy + age + precamp + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.force.minors.ses  <- cl(comp, faith.force.minors, comp$camp)       
      
  # Force - Birthplace
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "force", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.force.birthplace <- lm(govhelp ~ force*exposure + birthplace + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.force.birthplace.ses  <- cl(comp, faith.force.birthplace, comp$camp)
      
  # Force - LA
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "force", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      comp <- comp[comp$birthplace == "933",]
      
    # Model  
      faith.force.la <- lm(govhelp ~ force*exposure + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.force.la.ses  <- cl(comp, faith.force.la, comp$camp)    

 # Strikes
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.strikes <- lm(govhelp ~ strikes*exposure + precamp + age + gender + generation, data = comp)
    
    # Clustered Standard Errors 
      faith.strikes.ses  <- cl(comp, faith.strikes, comp$camp)
      
  # Strikes - Youth
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("govhelp", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.strikes.youth <- lm(govhelp ~ strikes*exposure + precamp + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.strikes.youth.ses  <- cl(comp, faith.strikes.youth, comp$camp)
      
  # Strikes - Age Interaction
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("govhelp", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.strikes.ageinter <- lm(govhelp ~ strikes*age + precamp  + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.strikes.ageinter.ses  <- cl(comp, faith.strikes.ageinter, comp$camp)    
      
  # Strikes - Minors
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("govhelp", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.strikes.minors <- lm(govhelp ~ strikes*minordummy + age + precamp  + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.strikes.minors.ses  <- cl(comp, faith.strikes.minors, comp$camp)
      
  # Strikes - Birthplace Controls
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "strikes", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.strikes.birthplace <- lm(govhelp ~ strikes*exposure + birthplace + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.strikes.birthplace.ses  <- cl(comp, faith.strikes.birthplace, comp$camp)  
      
  # Strikes - LA
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "strikes", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      comp <- comp[comp$birthplace == "933",]
      
    # Model  
      faith.strikes.la <- lm(govhelp ~ strikes*exposure + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.strikes.la.ses  <- cl(comp, faith.strikes.la, comp$camp)      
      
  # Militarism
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
    
    # Model  
      faith.militarism <- lm(govhelp ~ militarism*exposure + precamp + age + gender + generation, data = comp)
    
    # Clustered Standard Errors 
      faith.militarism.ses  <- cl(comp, faith.militarism, comp$camp)
      
  # Militarism Youth
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("govhelp", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.militarism.youth <- lm(govhelp ~ militarism*exposure + precamp + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.militarism.youth.ses  <- cl(comp, faith.militarism.youth, comp$camp)
      
  # Militarism Age Interaction
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("govhelp", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.militarism.ageinter <- lm(govhelp ~ militarism*age + precamp + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.militarism.ageinter.ses  <- cl(comp, faith.militarism.ageinter, comp$camp)    
      
  # Militarism Minors
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("govhelp", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.militarism.minors <- lm(govhelp ~ militarism*minordummy + age + precamp + gender + generation, data = comp)
      
      # Clustered Standard Errors 
      faith.militarism.minors.ses  <- cl(comp, faith.militarism.minors, comp$camp)   
      
  # Militarism - Birthplace
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "militarism", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      faith.militarism.birthplace <- lm(govhelp ~ militarism*exposure + birthplace + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.militarism.birthplace.ses  <- cl(comp, faith.militarism.birthplace, comp$camp)
      
  # Militarism - LA
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("govhelp", "militarism", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      comp <- comp[comp$birthplace == "933",]
      
    # Model  
      faith.militarism.la <- lm(govhelp ~ militarism*exposure + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      faith.militarism.la.ses  <- cl(comp, faith.militarism.la, comp$camp)    
      
################################################################################################################################################################################################################################################
# Thoughts on the Correct Type of Leadership
################################################################################################################################################################################################################################################
# Generate Combined Leadership Questions
  full$leader <- ifelse(full$generation == "Nisei", full$n.leader, ifelse(full$generation == "Sansei", full$s.leader, NA))
  full$leader <- ifelse(full$leader == "Protest", -1, ifelse(full$leader == "Orderly and comfortable", 1, ifelse(full$leader == "Can't generalize, both, neither" | full$leader == "Can't generalize; both; neither", 0, NA)))

# Run Regressions
  # Violence
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.violence <- lm(leader ~ violence*exposure + precamp + gender + age + generation, data = comp)  
    
    # Clustered Standard Errors 
      leader.violence.ses  <- cl(comp, leader.violence, comp$camp)
      
  # Violence Youth
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("leader", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.violence.youth <- lm(leader ~ violence*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.violence.youth.ses  <- cl(comp, leader.violence.youth, comp$camp)
      
  # Violence - Age Interaction
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("leader", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.violence.ageinter <- lm(leader ~ violence*age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.violence.ageinter.ses  <- cl(comp, leader.violence.ageinter, comp$camp)  
      
  # Violence - Minors
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("leader", "violence", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.violence.minors <- lm(leader ~ violence*minordummy + age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.violence.minors.ses  <- cl(comp, leader.violence.minors, comp$camp)       
      
  # Violence - Birthplace Controls
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "violence", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.violence.birthplace <- lm(leader ~ violence*exposure + birthplace + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.violence.birthplace.ses  <- cl(comp, leader.violence.birthplace, comp$camp)
      
  # Violence - LA
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "violence", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      comp <- comp[comp$birthplace == "933",]
      
    # Model  
      leader.violence.la <- lm(leader ~ violence*exposure + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.violence.la.ses  <- cl(comp, leader.violence.la, comp$camp)      
      
  # Force
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.force <- lm(leader ~ force*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.force.ses  <- cl(comp, leader.force, comp$camp)
      
  # Force - Youth
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("leader", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.force.youth <- lm(leader ~ force*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.force.youth.ses  <- cl(comp, leader.force.youth, comp$camp)
      
  # Force - Age Interaction
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("leader", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.force.ageinter <- lm(leader ~ force*age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.force.ageinter.ses  <- cl(comp, leader.force.ageinter, comp$camp)    
      
  # Force - Minors
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("leader", "force", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.force.minors <- lm(leader ~ force*minordummy + age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.force.minors.ses  <- cl(comp, leader.force.minors, comp$camp)   
      
  # Force - Birthplace Controls
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "force", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.force.birthplace <- lm(leader ~ force*exposure + birthplace + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.force.birthplace.ses  <- cl(comp, leader.force.birthplace, comp$camp)    
      
  # Force - LA Only
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "force", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      comp <- comp[comp$birthplace == "933",]
      
    # Model  
      leader.force.la <- lm(leader ~ force*exposure + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.force.la.ses  <- cl(comp, leader.force.la, comp$camp)    
      
  # Strikes 
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.strikes <- lm(leader ~ strikes*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.strikes.ses  <- cl(comp, leader.strikes, comp$camp)
      
  # Strikes - Youth
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("leader", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.strikes.youth <- lm(leader ~ strikes*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.strikes.youth.ses  <- cl(comp, leader.strikes.youth, comp$camp)
      
  # Strikes - Age Interaction
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("leader", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.strikes.ageinter <- lm(leader ~ strikes*age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.strikes.ageinter.ses  <- cl(comp, leader.strikes.ageinter, comp$camp)  
      
  # Strikes - Minors
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("leader", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.strikes.minors <- lm(leader ~ strikes*minordummy + age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.strikes.minors.ses  <- cl(comp, leader.strikes.minors, comp$camp)
      
  # Strikes - Birthplace
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "strikes", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.strikes.birthplace <- lm(leader ~ strikes*exposure + birthplace + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.strikes.birthplace.ses  <- cl(comp, leader.strikes.birthplace, comp$camp)  
      
  # Strikes - LA Only
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "strikes", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      comp <- comp[comp$birthplace == "933",]
      
    # Model  
      leader.strikes.la <- lm(leader ~ strikes*exposure + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.strikes.la.ses  <- cl(comp, leader.strikes.la, comp$camp)      
      
  # Militarism
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.militarism <- lm(leader ~ militarism*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.militarism.ses  <- cl(comp, leader.militarism, comp$camp)
      
  # Militarism - Youth
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("leader", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.militarism.youth <- lm(leader ~ militarism*exposure + precamp + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.militarism.youth.ses  <- cl(comp, leader.militarism.youth, comp$camp)    
      
  # Militarism - Age Interaction
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("leader", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.militarism.ageinter <- lm(leader ~ militarism*age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.militarism.ageinter.ses  <- cl(comp, leader.militarism.ageinter, comp$camp) 
      
  # Militarism - Minors
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned"), c("leader", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.militarism.minors <- lm(leader ~ militarism*minordummy + age + precamp + gender + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.militarism.minors.ses  <- cl(comp, leader.militarism.minors, comp$camp)     
      
  # Militarism - Birthplace Controls
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "militarism", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      leader.militarism.birthplace <- lm(leader ~ militarism*exposure + birthplace + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.militarism.birthplace.ses  <- cl(comp, leader.militarism.birthplace, comp$camp)    
      
  # Militarism - LA Only
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("leader", "militarism", "birthplace", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      comp <- comp[comp$birthplace == "933",]
      
    # Model  
      leader.militarism.la <- lm(leader ~ militarism*exposure + gender + age + generation, data = comp)  
      
    # Clustered Standard Errors 
      leader.militarism.la.ses  <- cl(comp, leader.militarism.la, comp$camp)         
      
################################################################################################################################################################################################################################################
# Being Approached for Political Advice
################################################################################################################################################################################################################################################
# Only asked of Sansei and Nisei
  full$poladvc <- ifelse(full$generation == "Nisei", full$n.poladv,
                       ifelse(full$generation == "Sansei", full$s.poladv, NA))
  
  full$poladvc <- ifelse(full$poladvc == "Yes", 1,
                       ifelse(full$poladvc == "No", 0, NA))  

# Violence
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.violence <- lm(poladvc ~ violence*exposure + precamp + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.violence.ses  <- cl(comp, polad.violence, comp$camp)
    
# Violence - Youth
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("poladvc", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.violence.youth <- lm(poladvc ~ violence*exposure + precamp + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.violence.youth.ses  <- cl(comp, polad.violence.youth, comp$camp)
    
# Violence Age Interaction
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned"), c("poladvc", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.violence.ageinter <- lm(poladvc ~ violence*age + precamp + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.violence.ageinter.ses  <- cl(comp, polad.violence.ageinter, comp$camp)    
    
# Violence Minors
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned"), c("poladvc", "violence", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.violence.minors <- lm(poladvc ~ violence*minordummy + age + precamp + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.violence.minors.ses  <- cl(comp, polad.violence.minors, comp$camp)     
    
# Violence - Birthplace
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "violence", "birthplace", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.violence.birthplace <- lm(poladvc ~ violence*exposure + birthplace + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.violence.birthplace.ses  <- cl(comp, polad.violence.birthplace, comp$camp)
    
# Violence - LA
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "violence", "birthplace", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    comp <- comp[comp$birthplace == "933",]
    
  # Model  
    polad.violence.la <- lm(poladvc ~ violence*exposure + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.violence.la.ses  <- cl(comp, polad.violence.la, comp$camp)    
    
# Force
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
  
  # Model  
    polad.force <- lm(poladvc ~ force*exposure + precamp + age + gender + generation, data = comp)
  
  # Clustered Standard Errors 
    polad.force.ses  <- cl(comp, polad.force, comp$camp)
    
# Force - Youth
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("poladvc", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.force.youth <- lm(poladvc ~ force*exposure + precamp + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.force.youth.ses  <- cl(comp, polad.force.youth, comp$camp) 
    
# Force - Age Interaction
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned"), c("poladvc", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.force.ageinter <- lm(poladvc ~ force*age + precamp + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.force.ageinter.ses  <- cl(comp, polad.force.ageinter, comp$camp)    
    
# Force - Minors
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned"), c("poladvc", "force", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.force.minors <- lm(poladvc ~ force*minordummy + age + precamp + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.force.minors.ses  <- cl(comp, polad.force.minors, comp$camp)  
    
# Force Birthplace
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "force", "birthplace", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.force.birthplace <- lm(poladvc ~ force*exposure + birthplace + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.force.birthplace.ses  <- cl(comp, polad.force.birthplace, comp$camp)
    
# Force Birthplace
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "force", "birthplace", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    comp <- comp[comp$birthplace == "933",]
    
  # Model  
    polad.force.la <- lm(poladvc ~ force*exposure + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.force.la.ses  <- cl(comp, polad.force.la, comp$camp)    
    
# Strikes
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.strikes <- lm(poladvc ~ strikes*exposure + precamp + age + gender + generation, data = comp)
  
  # Clustered Standard Errors 
    polad.strikes.ses  <- cl(comp, polad.strikes, comp$camp)
    
# Strikes - Youth
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("poladvc", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.strikes.youth <- lm(poladvc ~ strikes*exposure + precamp + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.strikes.youth.ses  <- cl(comp, polad.strikes.youth, comp$camp)
    

# Strikes Age Interaction
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned"), c("poladvc", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.strikes.ageinter <- lm(poladvc ~ strikes*age + precamp  + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.strikes.ageinter.ses  <- cl(comp, polad.strikes.ageinter, comp$camp)  
    
# Strikes Minors
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned"), c("poladvc", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.strikes.minors <- lm(poladvc ~ strikes*minordummy + age + precamp  + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.strikes.minors.ses  <- cl(comp, polad.strikes.minors, comp$camp)
    
# Strikes Birthplace
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "strikes", "birthplace", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.strikes.birthplace <- lm(poladvc ~ strikes*exposure + birthplace + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.strikes.birthplace.ses  <- cl(comp, polad.strikes.birthplace, comp$camp)
    
# Strikes LA
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "strikes", "birthplace", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    comp <- comp[comp$birthplace == "933",]
    
  # Model  
    polad.strikes.la <- lm(poladvc ~ strikes*exposure + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.strikes.la.ses  <- cl(comp, polad.strikes.la, comp$camp)    
    
# Militarism
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.militarism <- lm(poladvc ~ militarism*exposure + precamp + age + gender + generation, data = comp)
  
  # Clustered Standard Errors 
    polad.militarism.ses  <- cl(comp, polad.militarism, comp$camp)
    
# Militarism - Youth
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only") & full$age - 22 < 17, c("poladvc", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.militarism.youth <- lm(poladvc ~ militarism*exposure + precamp + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.militarism.youth.ses  <- cl(comp, polad.militarism.youth, comp$camp)  
    
# Militarism - Age Interaction
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned"), c("poladvc", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.militarism.ageinter <- lm(poladvc ~ militarism*age+ precamp  + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.militarism.ageinter.ses  <- cl(comp, polad.militarism.ageinter, comp$camp) 
    
# Militarism - Minors
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned"), c("poladvc", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure", "minordummy")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.militarism.minors <- lm(poladvc ~ militarism*minordummy + age+ precamp  + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.militarism.minors.ses  <- cl(comp, polad.militarism.minors, comp$camp)
    
# Militarism Birthplace
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "militarism", "birthplace", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    polad.militarism.birthplace <- lm(poladvc ~ militarism*exposure + birthplace + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.militarism.birthplace.ses  <- cl(comp, polad.militarism.birthplace, comp$camp)
    
# Militarism LA
  # Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("poladvc", "militarism", "birthplace", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    comp <- comp[comp$birthplace == "933",]
    
  # Model  
    polad.militarism.la <- lm(poladvc ~ militarism*exposure + age + gender + generation, data = comp)
    
  # Clustered Standard Errors 
    polad.militarism.la.ses  <- cl(comp, polad.militarism.la, comp$camp)    

################################################################################################################################################################################################################################################
# Party Affiliation and Vote Choice
################################################################################################################################################################################################################################################
# Code Vote Choice Across Generations (Note that I'm N/Aing out people who aren't old enough to vote or aren't citizens)    
  parties <- data.frame(party = unique(c(full$i.partyreg, full$n.party, full$s.party)), party.code = NA)
  parties$party.code <- ifelse(parties$party %in% c("Democrat", "Democratic"), -1,
                               ifelse(parties$party %in% c("Republican"), 1,
                                      ifelse(parties$party %in% c("Other", "Independant", "Independent", "Other", "Other (specify)"), 0, NA)))
  
    
  full$party <- ifelse(full$generation == "Issei", NA,
                       ifelse(full$generation == "Nisei", full$n.party,
                              ifelse(full$generation == "Sansei", full$s.party, NA)))
  
    
  full <- merge(full, parties, by = "party")
    
  
# Run Models  
  #Violence  
    #Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("party.code", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
    
    # Model  
      votechoice.violence <- lm(party.code ~ violence*exposure + precamp + age + gender + generation, data = comp)
  
    # Clustered Standard Errors 
      votechoice.violence.ses  <- cl(comp, votechoice.violence, comp$camp)
  
  # Force  
    # Complete Data for Standard Errors
      comp <-  full[full$exposure %in% c("self.interned", "fam.only"), c("party.code", "force", "precamp", "age", "gender", "camp", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      votechoice.force <- lm(party.code ~ force*exposure + precamp + age + gender, data = comp)
      
    # Clustered Standard Errors 
      votechoice.force.ses  <- cl(comp, votechoice.force, comp$camp)
      
  # Strikes  
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("party.code", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      votechoice.strikes <- lm(party.code ~ strikes*exposure + precamp + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      votechoice.strikes.ses  <- cl(comp, votechoice.strikes, comp$camp)
  
  # Militarism  
    # Complete Data for Standard Errors
      comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("party.code", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
      comp <- comp[complete.cases(comp),]
      
    # Model  
      votechoice.militarism <- lm(party.code ~ militarism*exposure + precamp + age + gender + generation, data = comp)
      
    # Clustered Standard Errors 
      votechoice.militarism.ses  <- cl(comp, votechoice.militarism, comp$camp)
      
      
      party.res <- data.frame(results = c(votechoice.violence$coefficients["violence"] + votechoice.violence$coefficients["violence:exposureself.interned"], votechoice.violence$coefficients["violence"],
                                          votechoice.strikes$coefficients["strikes"] + votechoice.strikes$coefficients["strikes:exposureself.interned"], votechoice.strikes$coefficients["strikes"],
                                          votechoice.force$coefficients["force"] + votechoice.force$coefficients["force:exposureself.interned"], votechoice.force$coefficients["force"],
                                          votechoice.militarism$coefficients["militarism"] + votechoice.militarism$coefficients["militarism:exposureself.interned"], votechoice.militarism$coefficients["militarism"]),
                                 ses = c(sqrt(votechoice.violence.ses[[2]]["violence", "violence"] + votechoice.violence.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * votechoice.violence.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                         sqrt(votechoice.violence.ses[[2]]["violence", "violence"]),
                                         sqrt(votechoice.strikes.ses[[2]]["strikes", "strikes"] + votechoice.strikes.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * votechoice.strikes.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                         sqrt(votechoice.strikes.ses[[2]]["strikes", "strikes"]),
                                         sqrt(votechoice.force.ses[[2]]["force", "force"] + votechoice.force.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * votechoice.force.ses[[2]]["force", "force:exposureself.interned"]) , 
                                         sqrt(votechoice.force.ses[[2]]["force", "force"]),
                                         sqrt(votechoice.militarism.ses[[2]]["militarism", "militarism"] + votechoice.militarism.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * votechoice.militarism.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                         sqrt(votechoice.militarism.ses[[2]]["militarism", "militarism"])))
      
      party.res$ub <- party.res$results + qnorm(.975) * party.res$ses
      party.res$lb <- party.res$results - qnorm(.975) * party.res$ses
      party.res$labels <- c(rep("Violence", 2), rep("Strikes", 2), rep("Force", 2), rep("Militarism", 2))
      party.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
      party.res$labels <- factor(party.res$labels, levels = c("Violence", "Strikes", "Force", "Militarism"), ordered = T)
      party.res$groups <- factor(party.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
      
    # Appendix J: Figure 22  
      pdf("si_fig22.pdf")
      ggplot(data = party.res, aes(x = labels, y = results, color = groups)) + 
        geom_hline(yintercept = 0, alpha = .5) +
        geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
        geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
        theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) + 
        labs(x = "Camp Characteristic", y = "More Likely to Support Democratic Party (-1) or Republican Party (1)", title = "Effects of Internment Camp Experience on Partisanship", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
      dev.off()
      
  # Voting is only Asked of Issei. Worth Noting that Most of them Vote Conditional On Being Citizens
      table(full$i.voted6064) 
      
  # Citizenship as an outcome
      full$citizen.num <- ifelse(full$i.citizen %in% c("1945", "1946", "1947", "1948", "1950", "1951", "1952", "1953", "1954", "1955", "1956", "1957", "1958", "1959", "1960", "1961", "Yes, no date given"), 1, 
                                 ifelse(full$i.citizen %in% c("DNA; no citizenship"), 0, NA))
                                 
    #Violence  
      #Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("citizen.num", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
      
      # Model  
        citizen.violence <- lm(citizen.num ~ violence + precamp + age + gender, data = comp)
      
      # Clustered Standard Errors 
        citizen.violence.ses  <- cl(comp, citizen.violence, comp$camp)
    
      #Force
        #Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("citizen.num", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
        # Model  
        citizen.force <- lm(citizen.num ~ force + precamp + age + gender, data = comp)
        
        # Clustered Standard Errors 
        citizen.force.ses  <- cl(comp, citizen.force, comp$camp)
      
        
      #Strikes
        #Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("citizen.num", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
        # Model  
        citizen.strikes <- lm(citizen.num ~ strikes + precamp + age + gender, data = comp)
        
        # Clustered Standard Errors 
        citizen.strikes.ses  <- cl(comp, citizen.strikes, comp$camp)   
        
      #Militarism
        #Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("citizen.num", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
        # Model  
        citizen.militarism <- lm(citizen.num ~ militarism + precamp + age + gender, data = comp)
        
        # Clustered Standard Errors 
        citizen.militarism.ses  <- cl(comp, citizen.militarism, comp$camp)   
        
        
        full$interned.numeric.num <- ifelse(full$interned.numeric == "Yes", 1, 
                                          ifelse(full$interned.numeric == "No", 0, NA))
      
        full$interned.numeric.num <- ifelse(full$exposure == "self.interned", 1, 0)
        
        
        citizen.res <- data.frame(results = c(citizen.violence$coefficients["violence"],
                                            citizen.strikes$coefficients["strikes"],
                                            citizen.force$coefficients["force"],
                                            citizen.militarism$coefficients["militarism"]),
                                ses = c(sqrt(citizen.violence.ses[[2]]["violence", "violence"]),
                                        sqrt(citizen.strikes.ses[[2]]["strikes", "strikes"]),
                                        sqrt(citizen.force.ses[[2]]["force", "force"]),
                                        sqrt(citizen.militarism.ses[[2]]["militarism", "militarism"])))
        
        citizen.res$ub <- citizen.res$results + qnorm(.975) * citizen.res$ses
        citizen.res$lb <- citizen.res$results - qnorm(.975) * citizen.res$ses
        citizen.res$labels <- c(rep("Violence", 1), rep("Strikes", 1), rep("Force", 1), rep("Militarism", 1))
        citizen.res$labels <- factor(citizen.res$labels, levels = c("Violence", "Strikes", "Force", "Militarism"), ordered = T)
        
        
        pdf("si_fig21.pdf")
        ggplot(data = citizen.res, aes(x = labels, y = results)) + 
          geom_hline(yintercept = 0, alpha = .5) +
          geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
          geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
          theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) + 
          labs(x = "Camp Characteristic", y = "Likelihood of Obtaining Citizenship After 1941", title = "Effects of Internment Camp Experience on Obtaining Citizenship", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
        dev.off()
      
      #Interned at All  
        comp <- full[, c("citizen.num", "interned.numeric.num", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
        # Model  
        citizen.interned <- lm(citizen.num ~ interned.numeric.num + age + gender + precamp, data = comp)
        
      
################################################################################################################################################################################################################################################
# Japanese American Continuum
################################################################################################################################################################################################################################################
# Get distribution by camp
  maincamps <- full[full$camp %in%c("Rivers, AZ", "Amache, CO", "Heart Mountain, WY", "Denson, AR", "Manzanar, CA", "Hunt, ID", "Poston, AZ", "McGehee, AR","Topaz, UT", "Newell, CA"),]  
  
  pdf("si_fig23.pdf")    
      ggplot(maincamps, aes(x = maincamps$i.continuum, group = maincamps$campname)) +
      geom_bar() + facet_wrap(~campname) +
      scale_x_discrete(limit = c("Still fully Japanese (little or no change)", "About 80% Japanese",
                                 "About 65% Japanese", "About 50-50", "About 65% Americanized", "About 80% Americanized",
                                 "About 100% Americanized"),
                         labels = c("100% J","80% J","65% J", "50%-50%", "65% A", "80% A", "100% A")) + 
      theme(axis.text.x = element_text(angle = 60, hjust = 1)) + xlab("Position on the Japanese-American Continuum") + ylab("Frequency")
  dev.off()
  
  full$i.continuum.num <- ifelse(full$i.continuum == "Still fully Japanese (little or no change)", 1,
                                 ifelse(full$i.continuum == "About 80% Japanese", 2,
                                        ifelse(full$i.continuum == "About 65% Japanese", 3,
                                               ifelse(full$i.continuum == "About 50-50", 4,
                                                      ifelse(full$i.continuum == "About 65% Americanized", 5,
                                                             ifelse(full$i.continuum == "About 80% Americanized", 6, 
                                                                    ifelse(full$i.continuum == "About 100% Americanized", 7, NA)))))))
  
  
#Violence  
  #Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("i.continuum.num", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    continuum.violence <- lm(i.continuum.num ~ violence + precamp + age + gender, data = comp)
  
  # Clustered Standard Errors 
    continuum.violence.ses  <- cl(comp, continuum.violence, comp$camp)
  
#Force
  #Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("i.continuum.num", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    continuum.force <- lm(i.continuum.num ~ force + precamp + age + gender, data = comp)
  
  # Clustered Standard Errors 
    continuum.force.ses  <- cl(comp, continuum.force, comp$camp)
  
#Strikes
  #Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("i.continuum.num", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    continuum.strikes <- lm(i.continuum.num ~ strikes + precamp + age + gender, data = comp)
  
  # Clustered Standard Errors 
    continuum.strikes.ses  <- cl(comp, continuum.strikes, comp$camp)   
  
#Militarism
  #Complete Data for Standard Errors
    comp <- full[full$exposure %in% c("self.interned", "fam.only"), c("i.continuum.num", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
    comp <- comp[complete.cases(comp),]
    
  # Model  
    continuum.militarism <- lm(i.continuum.num ~ militarism + precamp + age + gender, data = comp)
  
  # Clustered Standard Errors 
    continuum.militarism.ses  <- cl(comp, continuum.militarism, comp$camp)   
  
  
    full$interned.numeric.num <- ifelse(full$interned.numeric == "Yes", 1, 
                                        ifelse(full$interned.numeric == "No", 0, NA))
    
    full$interned.numeric.num <- ifelse(full$exposure == "self.interned", 1, 0)
    
  
    continuum.res <- data.frame(results = c(continuum.violence$coefficients["violence"],
                                          continuum.strikes$coefficients["strikes"],
                                          continuum.force$coefficients["force"],
                                          continuum.militarism$coefficients["militarism"]),
                              ses = c(sqrt(continuum.violence.ses[[2]]["violence", "violence"]),
                                      sqrt(continuum.strikes.ses[[2]]["strikes", "strikes"]),
                                      sqrt(continuum.force.ses[[2]]["force", "force"]),
                                      sqrt(continuum.militarism.ses[[2]]["militarism", "militarism"])))
    
    continuum.res$ub <- continuum.res$results + qnorm(.975) * continuum.res$ses
    continuum.res$lb <- continuum.res$results - qnorm(.975) * continuum.res$ses
    continuum.res$labels <- c(rep("Violence", 1), rep("Strikes", 1), rep("Force", 1), rep("Militarism", 1))
    continuum.res$labels <- factor(continuum.res$labels, levels = c("Violence", "Strikes", "Force", "Militarism"), ordered = T)
    
    pdf("si_fig24.pdf")
    ggplot(data = continuum.res, aes(x = labels, y = results)) + 
      geom_hline(yintercept = 0, alpha = .5) +
      geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
      geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
      theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) + 
      labs(x = "Camp Characteristic", y = "Position on Japanese (Left) - American (Right) Continuum", title = "Effects of Internment Camp Experience on Assimilation", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
    dev.off()
  
################################################################################################################################################################################################################################################
# Plots by Treatment
################################################################################################################################################################################################################################################
# Effects of Violence
  violence.res <- data.frame(results = c(interest.violence$coefficients["violence"] + interest.violence$coefficients["violence:exposureself.interned"], interest.violence$coefficients["violence"],
                                         faith.violence$coefficients["violence"] + faith.violence$coefficients["violence:exposureself.interned"], faith.violence$coefficients["violence"],
                                         leader.violence$coefficients["violence"] + leader.violence$coefficients["violence:exposureself.interned"], leader.violence$coefficients["violence"], 
                                         polad.violence$coefficients["violence"] + polad.violence$coefficients["violence:exposureself.interned"], polad.violence$coefficients["violence"]),
                          ses = c(sqrt(interest.violence.ses[[2]]["violence", "violence"] + interest.violence.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * interest.violence.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                  sqrt(interest.violence.ses[[2]]["violence", "violence"]),
                                  sqrt(faith.violence.ses[[2]]["violence", "violence"] + faith.violence.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * faith.violence.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                  sqrt(faith.violence.ses[[2]]["violence", "violence"]), 
                                  sqrt(leader.violence.ses[[2]]["violence", "violence"] + leader.violence.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * leader.violence.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                  sqrt(leader.violence.ses[[2]]["violence", "violence"]),
                                  sqrt(polad.violence.ses[[2]]["violence", "violence"] + polad.violence.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * polad.violence.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                  sqrt(polad.violence.ses[[2]]["violence", "violence"])))
  violence.res$ub <- violence.res$results + qnorm(.975) * violence.res$ses
  violence.res$lb <- violence.res$results - qnorm(.975) * violence.res$ses
  violence.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  violence.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  violence.res$labels <- factor(violence.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  violence.res$groups <- factor(violence.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  violence.youth.res <- data.frame(results = c(interest.violence.youth$coefficients["violence"] + interest.violence.youth$coefficients["violence:exposureself.interned"], interest.violence.youth$coefficients["violence"],
                                         faith.violence.youth$coefficients["violence"] + faith.violence.youth$coefficients["violence:exposureself.interned"], faith.violence.youth$coefficients["violence"],
                                         leader.violence.youth$coefficients["violence"] + leader.violence.youth$coefficients["violence:exposureself.interned"], leader.violence.youth$coefficients["violence"], 
                                         polad.violence.youth$coefficients["violence"] + polad.violence.youth$coefficients["violence:exposureself.interned"], polad.violence.youth$coefficients["violence"]),
                             ses = c(sqrt(interest.violence.youth.ses[[2]]["violence", "violence"] + interest.violence.youth.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * interest.violence.youth.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(interest.violence.youth.ses[[2]]["violence", "violence"]),
                                     sqrt(faith.violence.youth.ses[[2]]["violence", "violence"] + faith.violence.youth.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * faith.violence.youth.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(faith.violence.youth.ses[[2]]["violence", "violence"]), 
                                     sqrt(leader.violence.youth.ses[[2]]["violence", "violence"] + leader.violence.youth.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * leader.violence.youth.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(leader.violence.youth.ses[[2]]["violence", "violence"]),
                                     sqrt(polad.violence.youth.ses[[2]]["violence", "violence"] + polad.violence.youth.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * polad.violence.youth.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(polad.violence.youth.ses[[2]]["violence", "violence"])))
  violence.youth.res$ub <- violence.youth.res$results + qnorm(.975) * violence.youth.res$ses
  violence.youth.res$lb <- violence.youth.res$results - qnorm(.975) * violence.youth.res$ses
  violence.youth.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  violence.youth.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  violence.youth.res$labels <- factor(violence.youth.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  violence.youth.res$groups <- factor(violence.youth.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  
  violence.la.res <- data.frame(results = c(interest.violence.la$coefficients["violence"] + interest.violence.la$coefficients["violence:exposureself.interned"], interest.violence.la$coefficients["violence"],
                                         faith.violence.la$coefficients["violence"] + faith.violence.la$coefficients["violence:exposureself.interned"], faith.violence.la$coefficients["violence"],
                                         leader.violence.la$coefficients["violence"] + leader.violence.la$coefficients["violence:exposureself.interned"], leader.violence.la$coefficients["violence"], 
                                         polad.violence.la$coefficients["violence"] + polad.violence.la$coefficients["violence:exposureself.interned"], polad.violence.la$coefficients["violence"]),
                             ses = c(sqrt(interest.violence.la.ses[[2]]["violence", "violence"] + interest.violence.la.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * interest.violence.la.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(interest.violence.la.ses[[2]]["violence", "violence"]),
                                     sqrt(faith.violence.la.ses[[2]]["violence", "violence"] + faith.violence.la.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * faith.violence.la.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(faith.violence.la.ses[[2]]["violence", "violence"]), 
                                     sqrt(leader.violence.la.ses[[2]]["violence", "violence"] + leader.violence.la.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * leader.violence.la.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(leader.violence.la.ses[[2]]["violence", "violence"]),
                                     sqrt(polad.violence.la.ses[[2]]["violence", "violence"] + polad.violence.la.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * polad.violence.la.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(polad.violence.la.ses[[2]]["violence", "violence"])))
  violence.la.res$ub <- violence.la.res$results + qnorm(.975) * violence.la.res$ses
  violence.la.res$lb <- violence.la.res$results - qnorm(.975) * violence.la.res$ses
  violence.la.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  violence.la.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  violence.la.res$labels <- factor(violence.la.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  violence.la.res$groups <- factor(violence.la.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  violence.birthplace.res <- data.frame(results = c(interest.violence.birthplace$coefficients["violence"] + interest.violence.birthplace$coefficients["violence:exposureself.interned"], interest.violence.birthplace$coefficients["violence"],
                                         faith.violence.birthplace$coefficients["violence"] + faith.violence.birthplace$coefficients["violence:exposureself.interned"], faith.violence.birthplace$coefficients["violence"],
                                         leader.violence.birthplace$coefficients["violence"] + leader.violence.birthplace$coefficients["violence:exposureself.interned"], leader.violence.birthplace$coefficients["violence"], 
                                         polad.violence.birthplace$coefficients["violence"] + polad.violence.birthplace$coefficients["violence:exposureself.interned"], polad.violence.birthplace$coefficients["violence"]),
                             ses = c(sqrt(interest.violence.birthplace.ses[[2]]["violence", "violence"] + interest.violence.birthplace.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * interest.violence.birthplace.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(interest.violence.birthplace.ses[[2]]["violence", "violence"]),
                                     sqrt(faith.violence.birthplace.ses[[2]]["violence", "violence"] + faith.violence.birthplace.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * faith.violence.birthplace.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(faith.violence.birthplace.ses[[2]]["violence", "violence"]), 
                                     sqrt(leader.violence.birthplace.ses[[2]]["violence", "violence"] + leader.violence.birthplace.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * leader.violence.birthplace.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(leader.violence.birthplace.ses[[2]]["violence", "violence"]),
                                     sqrt(polad.violence.birthplace.ses[[2]]["violence", "violence"] + polad.violence.birthplace.ses[[2]]["violence:exposureself.interned", "violence:exposureself.interned"] + 2 * polad.violence.birthplace.ses[[2]]["violence", "violence:exposureself.interned"]) , 
                                     sqrt(polad.violence.birthplace.ses[[2]]["violence", "violence"])))
  violence.birthplace.res$ub <- violence.birthplace.res$results + qnorm(.975) * violence.birthplace.res$ses
  violence.birthplace.res$lb <- violence.birthplace.res$results - qnorm(.975) * violence.birthplace.res$ses
  violence.birthplace.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  violence.birthplace.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  violence.birthplace.res$labels <- factor(violence.birthplace.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  violence.birthplace.res$groups <- factor(violence.birthplace.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  
  
  pdf("ms_fig7.pdf")
  ggplot(data = violence.res, aes(x = labels, y = results, color = groups)) + 
    geom_hline(yintercept = 0, alpha = .5) + 
    geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
    geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
    theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) + 
    labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Violence While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig16b.pdf")
  ggplot(data = violence.youth.res, aes(x = labels, y = results, color = groups)) + 
    geom_hline(yintercept = 0, alpha = .5) + 
    geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
    geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
    theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) + 
    labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Violence While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig15b.pdf")
  ggplot(data = violence.la.res, aes(x = labels, y = results, color = groups)) + 
    geom_hline(yintercept = 0, alpha = .5) + 
    geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
    geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
    theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) + 
    labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Violence While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig14b.pdf")
  ggplot(data = violence.birthplace.res, aes(x = labels, y = results, color = groups)) + 
    geom_hline(yintercept = 0, alpha = .5) + 
    geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
    geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
    theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) + 
    labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Violence While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  
  
# Effects of Strikes  
  strikes.res <- data.frame(results = c(interest.strikes$coefficients["strikes"] + interest.strikes$coefficients["strikes:exposureself.interned"], interest.strikes$coefficients["strikes"],
                                         faith.strikes$coefficients["strikes"] + faith.strikes$coefficients["strikes:exposureself.interned"], faith.strikes$coefficients["strikes"],
                                         leader.strikes$coefficients["strikes"] + leader.strikes$coefficients["strikes:exposureself.interned"], leader.strikes$coefficients["strikes"], 
                                         polad.strikes$coefficients["strikes"] + polad.strikes$coefficients["strikes:exposureself.interned"], polad.strikes$coefficients["strikes"]),
                             ses = c(sqrt(interest.strikes.ses[[2]]["strikes", "strikes"] + interest.strikes.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * interest.strikes.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                     sqrt(interest.strikes.ses[[2]]["strikes", "strikes"]),
                                     sqrt(faith.strikes.ses[[2]]["strikes", "strikes"] + faith.strikes.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * faith.strikes.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                     sqrt(faith.strikes.ses[[2]]["strikes", "strikes"]), 
                                     sqrt(leader.strikes.ses[[2]]["strikes", "strikes"] + leader.strikes.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * leader.strikes.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                     sqrt(leader.strikes.ses[[2]]["strikes", "strikes"]),
                                     sqrt(polad.strikes.ses[[2]]["strikes", "strikes"] + polad.strikes.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * polad.strikes.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                     sqrt(polad.strikes.ses[[2]]["strikes", "strikes"])))
  strikes.res$ub <- strikes.res$results + qnorm(.975) * strikes.res$ses
  strikes.res$lb <- strikes.res$results - qnorm(.975) * strikes.res$ses
  strikes.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  strikes.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  strikes.res$labels <- factor(strikes.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  strikes.res$groups <- factor(strikes.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  strikes.youth.res <- data.frame(results = c(interest.strikes.youth$coefficients["strikes"] + interest.strikes.youth$coefficients["strikes:exposureself.interned"], interest.strikes.youth$coefficients["strikes"],
                                        faith.strikes.youth$coefficients["strikes"] + faith.strikes.youth$coefficients["strikes:exposureself.interned"], faith.strikes.youth$coefficients["strikes"],
                                        leader.strikes.youth$coefficients["strikes"] + leader.strikes.youth$coefficients["strikes:exposureself.interned"], leader.strikes.youth$coefficients["strikes"], 
                                        polad.strikes.youth$coefficients["strikes"] + polad.strikes.youth$coefficients["strikes:exposureself.interned"], polad.strikes.youth$coefficients["strikes"]),
                            ses = c(sqrt(interest.strikes.youth.ses[[2]]["strikes", "strikes"] + interest.strikes.youth.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * interest.strikes.youth.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                    sqrt(interest.strikes.youth.ses[[2]]["strikes", "strikes"]),
                                    sqrt(faith.strikes.youth.ses[[2]]["strikes", "strikes"] + faith.strikes.youth.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * faith.strikes.youth.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                    sqrt(faith.strikes.youth.ses[[2]]["strikes", "strikes"]), 
                                    sqrt(leader.strikes.youth.ses[[2]]["strikes", "strikes"] + leader.strikes.youth.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * leader.strikes.youth.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                    sqrt(leader.strikes.youth.ses[[2]]["strikes", "strikes"]),
                                    sqrt(polad.strikes.youth.ses[[2]]["strikes", "strikes"] + polad.strikes.youth.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * polad.strikes.youth.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                    sqrt(polad.strikes.youth.ses[[2]]["strikes", "strikes"])))
  strikes.youth.res$ub <- strikes.youth.res$results + qnorm(.975) * strikes.youth.res$ses
  strikes.youth.res$lb <- strikes.youth.res$results - qnorm(.975) * strikes.youth.res$ses
  strikes.youth.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  strikes.youth.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  strikes.youth.res$labels <- factor(strikes.youth.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  strikes.youth.res$groups <- factor(strikes.youth.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  strikes.la.res <- data.frame(results = c(interest.strikes.la$coefficients["strikes"] + interest.strikes.la$coefficients["strikes:exposureself.interned"], interest.strikes.la$coefficients["strikes"],
                                              faith.strikes.la$coefficients["strikes"] + faith.strikes.la$coefficients["strikes:exposureself.interned"], faith.strikes.la$coefficients["strikes"],
                                              leader.strikes.la$coefficients["strikes"] + leader.strikes.la$coefficients["strikes:exposureself.interned"], leader.strikes.la$coefficients["strikes"], 
                                              polad.strikes.la$coefficients["strikes"] + polad.strikes.la$coefficients["strikes:exposureself.interned"], polad.strikes.la$coefficients["strikes"]),
                                  ses = c(sqrt(interest.strikes.la.ses[[2]]["strikes", "strikes"] + interest.strikes.la.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * interest.strikes.la.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                          sqrt(interest.strikes.la.ses[[2]]["strikes", "strikes"]),
                                          sqrt(faith.strikes.la.ses[[2]]["strikes", "strikes"] + faith.strikes.la.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * faith.strikes.la.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                          sqrt(faith.strikes.la.ses[[2]]["strikes", "strikes"]), 
                                          sqrt(leader.strikes.la.ses[[2]]["strikes", "strikes"] + leader.strikes.la.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * leader.strikes.la.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                          sqrt(leader.strikes.la.ses[[2]]["strikes", "strikes"]),
                                          sqrt(polad.strikes.la.ses[[2]]["strikes", "strikes"] + polad.strikes.la.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * polad.strikes.la.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                          sqrt(polad.strikes.la.ses[[2]]["strikes", "strikes"])))
  strikes.la.res$ub <- strikes.la.res$results + qnorm(.975) * strikes.la.res$ses
  strikes.la.res$lb <- strikes.la.res$results - qnorm(.975) * strikes.la.res$ses
  strikes.la.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  strikes.la.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  strikes.la.res$labels <- factor(strikes.la.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  strikes.la.res$groups <- factor(strikes.la.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  strikes.birthplace.res <- data.frame(results = c(interest.strikes.birthplace$coefficients["strikes"] + interest.strikes.birthplace$coefficients["strikes:exposureself.interned"], interest.strikes.birthplace$coefficients["strikes"],
                                        faith.strikes.birthplace$coefficients["strikes"] + faith.strikes.birthplace$coefficients["strikes:exposureself.interned"], faith.strikes.birthplace$coefficients["strikes"],
                                        leader.strikes.birthplace$coefficients["strikes"] + leader.strikes.birthplace$coefficients["strikes:exposureself.interned"], leader.strikes.birthplace$coefficients["strikes"], 
                                        polad.strikes.birthplace$coefficients["strikes"] + polad.strikes.birthplace$coefficients["strikes:exposureself.interned"], polad.strikes.birthplace$coefficients["strikes"]),
                            ses = c(sqrt(interest.strikes.birthplace.ses[[2]]["strikes", "strikes"] + interest.strikes.birthplace.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * interest.strikes.birthplace.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                    sqrt(interest.strikes.birthplace.ses[[2]]["strikes", "strikes"]),
                                    sqrt(faith.strikes.birthplace.ses[[2]]["strikes", "strikes"] + faith.strikes.birthplace.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * faith.strikes.birthplace.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                    sqrt(faith.strikes.birthplace.ses[[2]]["strikes", "strikes"]), 
                                    sqrt(leader.strikes.birthplace.ses[[2]]["strikes", "strikes"] + leader.strikes.birthplace.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * leader.strikes.birthplace.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                    sqrt(leader.strikes.birthplace.ses[[2]]["strikes", "strikes"]),
                                    sqrt(polad.strikes.birthplace.ses[[2]]["strikes", "strikes"] + polad.strikes.birthplace.ses[[2]]["strikes:exposureself.interned", "strikes:exposureself.interned"] + 2 * polad.strikes.birthplace.ses[[2]]["strikes", "strikes:exposureself.interned"]) , 
                                    sqrt(polad.strikes.birthplace.ses[[2]]["strikes", "strikes"])))
  strikes.birthplace.res$ub <- strikes.birthplace.res$results + qnorm(.975) * strikes.birthplace.res$ses
  strikes.birthplace.res$lb <- strikes.birthplace.res$results - qnorm(.975) * strikes.birthplace.res$ses
  strikes.birthplace.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  strikes.birthplace.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  strikes.birthplace.res$labels <- factor(strikes.birthplace.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  strikes.birthplace.res$groups <- factor(strikes.birthplace.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  
  pdf("ms_fig6.pdf")
  ggplot(data = strikes.res, aes(x = labels, y = results, color = groups)) + 
    geom_hline(yintercept = 0, alpha = .5) + 
    geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
    geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
    theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
    labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Demonstrations While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig16a.pdf")
  ggplot(data = strikes.youth.res, aes(x = labels, y = results, color = groups)) + 
    geom_hline(yintercept = 0, alpha = .5) + 
    geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
    geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
    theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
    labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Demonstrations While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig15a.pdf")
  ggplot(data = strikes.la.res, aes(x = labels, y = results, color = groups)) + 
    geom_hline(yintercept = 0, alpha = .5) + 
    geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
    geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
    theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
    labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Demonstrations While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig14a.pdf")
  ggplot(data = strikes.birthplace.res, aes(x = labels, y = results, color = groups)) + 
    geom_hline(yintercept = 0, alpha = .5) + 
    geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
    geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
    theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
    labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Demonstrations While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  
  
# Effects of Militarization
  militarism.res <- data.frame(results = c(interest.militarism$coefficients["militarism"] + interest.militarism$coefficients["militarism:exposureself.interned"], interest.militarism$coefficients["militarism"],
                                         faith.militarism$coefficients["militarism"] + faith.militarism$coefficients["militarism:exposureself.interned"], faith.militarism$coefficients["militarism"],
                                         leader.militarism$coefficients["militarism"] + leader.militarism$coefficients["militarism:exposureself.interned"], leader.militarism$coefficients["militarism"], 
                                         polad.militarism$coefficients["militarism"] + polad.militarism$coefficients["militarism:exposureself.interned"], polad.militarism$coefficients["militarism"]),
                             ses = c(sqrt(interest.militarism.ses[[2]]["militarism", "militarism"] + interest.militarism.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * interest.militarism.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                     sqrt(interest.militarism.ses[[2]]["militarism", "militarism"]),
                                     sqrt(faith.militarism.ses[[2]]["militarism", "militarism"] + faith.militarism.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * faith.militarism.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                     sqrt(faith.militarism.ses[[2]]["militarism", "militarism"]), 
                                     sqrt(leader.militarism.ses[[2]]["militarism", "militarism"] + leader.militarism.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * leader.militarism.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                     sqrt(leader.militarism.ses[[2]]["militarism", "militarism"]),
                                     sqrt(polad.militarism.ses[[2]]["militarism", "militarism"] + polad.militarism.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * polad.militarism.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                     sqrt(polad.militarism.ses[[2]]["militarism", "militarism"])))
  militarism.res$ub <- militarism.res$results + qnorm(.975) * militarism.res$ses
  militarism.res$lb <- militarism.res$results - qnorm(.975) * militarism.res$ses
  militarism.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  militarism.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  militarism.res$labels <- factor(militarism.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  militarism.res$groups <- factor(militarism.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  militarism.youth.res <- data.frame(results = c(interest.militarism.youth$coefficients["militarism"] + interest.militarism.youth$coefficients["militarism:exposureself.interned"], interest.militarism.youth$coefficients["militarism"],
                                           faith.militarism.youth$coefficients["militarism"] + faith.militarism.youth$coefficients["militarism:exposureself.interned"], faith.militarism.youth$coefficients["militarism"],
                                           leader.militarism.youth$coefficients["militarism"] + leader.militarism.youth$coefficients["militarism:exposureself.interned"], leader.militarism.youth$coefficients["militarism"], 
                                           polad.militarism.youth$coefficients["militarism"] + polad.militarism.youth$coefficients["militarism:exposureself.interned"], polad.militarism.youth$coefficients["militarism"]),
                               ses = c(sqrt(interest.militarism.youth.ses[[2]]["militarism", "militarism"] + interest.militarism.youth.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * interest.militarism.youth.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                       sqrt(interest.militarism.youth.ses[[2]]["militarism", "militarism"]),
                                       sqrt(faith.militarism.youth.ses[[2]]["militarism", "militarism"] + faith.militarism.youth.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * faith.militarism.youth.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                       sqrt(faith.militarism.youth.ses[[2]]["militarism", "militarism"]), 
                                       sqrt(leader.militarism.youth.ses[[2]]["militarism", "militarism"] + leader.militarism.youth.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * leader.militarism.youth.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                       sqrt(leader.militarism.youth.ses[[2]]["militarism", "militarism"]),
                                       sqrt(polad.militarism.youth.ses[[2]]["militarism", "militarism"] + polad.militarism.youth.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * polad.militarism.youth.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                       sqrt(polad.militarism.youth.ses[[2]]["militarism", "militarism"])))
  militarism.youth.res$ub <- militarism.youth.res$results + qnorm(.975) * militarism.youth.res$ses
  militarism.youth.res$lb <- militarism.youth.res$results - qnorm(.975) * militarism.youth.res$ses
  militarism.youth.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  militarism.youth.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  militarism.youth.res$labels <- factor(militarism.youth.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  militarism.youth.res$groups <- factor(militarism.youth.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  
  militarism.la.res <- data.frame(results = c(interest.militarism.la$coefficients["militarism"] + interest.militarism.la$coefficients["militarism:exposureself.interned"], interest.militarism.la$coefficients["militarism"],
                                                 faith.militarism.la$coefficients["militarism"] + faith.militarism.la$coefficients["militarism:exposureself.interned"], faith.militarism.la$coefficients["militarism"],
                                                 leader.militarism.la$coefficients["militarism"] + leader.militarism.la$coefficients["militarism:exposureself.interned"], leader.militarism.la$coefficients["militarism"], 
                                                 polad.militarism.la$coefficients["militarism"] + polad.militarism.la$coefficients["militarism:exposureself.interned"], polad.militarism.la$coefficients["militarism"]),
                                     ses = c(sqrt(interest.militarism.la.ses[[2]]["militarism", "militarism"] + interest.militarism.la.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * interest.militarism.la.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                             sqrt(interest.militarism.la.ses[[2]]["militarism", "militarism"]),
                                             sqrt(faith.militarism.la.ses[[2]]["militarism", "militarism"] + faith.militarism.la.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * faith.militarism.la.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                             sqrt(faith.militarism.la.ses[[2]]["militarism", "militarism"]), 
                                             sqrt(leader.militarism.la.ses[[2]]["militarism", "militarism"] + leader.militarism.la.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * leader.militarism.la.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                             sqrt(leader.militarism.la.ses[[2]]["militarism", "militarism"]),
                                             sqrt(polad.militarism.la.ses[[2]]["militarism", "militarism"] + polad.militarism.la.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * polad.militarism.la.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                             sqrt(polad.militarism.la.ses[[2]]["militarism", "militarism"])))
  militarism.la.res$ub <- militarism.la.res$results + qnorm(.975) * militarism.la.res$ses
  militarism.la.res$lb <- militarism.la.res$results - qnorm(.975) * militarism.la.res$ses
  militarism.la.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  militarism.la.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  militarism.la.res$labels <- factor(militarism.la.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  militarism.la.res$groups <- factor(militarism.la.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  militarism.birthplace.res <- data.frame(results = c(interest.militarism.birthplace$coefficients["militarism"] + interest.militarism.birthplace$coefficients["militarism:exposureself.interned"], interest.militarism.birthplace$coefficients["militarism"],
                                           faith.militarism.birthplace$coefficients["militarism"] + faith.militarism.birthplace$coefficients["militarism:exposureself.interned"], faith.militarism.birthplace$coefficients["militarism"],
                                           leader.militarism.birthplace$coefficients["militarism"] + leader.militarism.birthplace$coefficients["militarism:exposureself.interned"], leader.militarism.birthplace$coefficients["militarism"], 
                                           polad.militarism.birthplace$coefficients["militarism"] + polad.militarism.birthplace$coefficients["militarism:exposureself.interned"], polad.militarism.birthplace$coefficients["militarism"]),
                               ses = c(sqrt(interest.militarism.birthplace.ses[[2]]["militarism", "militarism"] + interest.militarism.birthplace.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * interest.militarism.birthplace.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                       sqrt(interest.militarism.birthplace.ses[[2]]["militarism", "militarism"]),
                                       sqrt(faith.militarism.birthplace.ses[[2]]["militarism", "militarism"] + faith.militarism.birthplace.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * faith.militarism.birthplace.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                       sqrt(faith.militarism.birthplace.ses[[2]]["militarism", "militarism"]), 
                                       sqrt(leader.militarism.birthplace.ses[[2]]["militarism", "militarism"] + leader.militarism.birthplace.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * leader.militarism.birthplace.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                       sqrt(leader.militarism.birthplace.ses[[2]]["militarism", "militarism"]),
                                       sqrt(polad.militarism.birthplace.ses[[2]]["militarism", "militarism"] + polad.militarism.birthplace.ses[[2]]["militarism:exposureself.interned", "militarism:exposureself.interned"] + 2 * polad.militarism.birthplace.ses[[2]]["militarism", "militarism:exposureself.interned"]) , 
                                       sqrt(polad.militarism.birthplace.ses[[2]]["militarism", "militarism"])))
  militarism.birthplace.res$ub <- militarism.birthplace.res$results + qnorm(.975) * militarism.birthplace.res$ses
  militarism.birthplace.res$lb <- militarism.birthplace.res$results - qnorm(.975) * militarism.birthplace.res$ses
  militarism.birthplace.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  militarism.birthplace.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  militarism.birthplace.res$labels <- factor(militarism.birthplace.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  militarism.birthplace.res$groups <- factor(militarism.birthplace.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  
  
  pdf("ms_fig9a.pdf")
  ggplot(data = militarism.res, aes(x = labels, y = results, color = groups)) + 
    geom_hline(yintercept = 0, alpha = .5) + 
    geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
    geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
    theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
    labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Being Interned in Highly Militarized Camp", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig16d.pdf")
    ggplot(data = militarism.youth.res, aes(x = labels, y = results, color = groups)) + 
      geom_hline(yintercept = 0, alpha = .5) + 
      geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
      geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
      theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
      labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Being Interned in Highly Militarized Camp", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig15d.pdf")
    ggplot(data = militarism.la.res, aes(x = labels, y = results, color = groups)) + 
      geom_hline(yintercept = 0, alpha = .5) + 
      geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
      geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
      theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
      labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Being Interned in Highly Militarized Camp", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig14d.pdf")
  ggplot(data = militarism.birthplace.res, aes(x = labels, y = results, color = groups)) + 
    geom_hline(yintercept = 0, alpha = .5) + 
    geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
    geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
    theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
    labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Being Interned in Highly Militarized Camp", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  
  
# Effects of Force
  force.res <- data.frame(results = c(interest.force$coefficients["force"] + interest.force$coefficients["force:exposureself.interned"], interest.force$coefficients["force"],
                                         faith.force$coefficients["force"] + faith.force$coefficients["force:exposureself.interned"], faith.force$coefficients["force"],
                                         leader.force$coefficients["force"] + leader.force$coefficients["force:exposureself.interned"], leader.force$coefficients["force"], 
                                         polad.force$coefficients["force"] + polad.force$coefficients["force:exposureself.interned"], polad.force$coefficients["force"]),
                             ses = c(sqrt(interest.force.ses[[2]]["force", "force"] + interest.force.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * interest.force.ses[[2]]["force", "force:exposureself.interned"]) , 
                                     sqrt(interest.force.ses[[2]]["force", "force"]),
                                     sqrt(faith.force.ses[[2]]["force", "force"] + faith.force.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * faith.force.ses[[2]]["force", "force:exposureself.interned"]) , 
                                     sqrt(faith.force.ses[[2]]["force", "force"]), 
                                     sqrt(leader.force.ses[[2]]["force", "force"] + leader.force.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * leader.force.ses[[2]]["force", "force:exposureself.interned"]) , 
                                     sqrt(leader.force.ses[[2]]["force", "force"]),
                                     sqrt(polad.force.ses[[2]]["force", "force"] + polad.force.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * polad.force.ses[[2]]["force", "force:exposureself.interned"]) , 
                                     sqrt(polad.force.ses[[2]]["force", "force"])))
  force.res$ub <- force.res$results + qnorm(.975) * force.res$ses
  force.res$lb <- force.res$results - qnorm(.975) * force.res$ses
  force.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  force.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  force.res$labels <- factor(force.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  force.res$groups <- factor(force.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  

  force.youth.res <- data.frame(results = c(interest.force.youth$coefficients["force"] + interest.force.youth$coefficients["force:exposureself.interned"], interest.force.youth$coefficients["force"],
                                      faith.force.youth$coefficients["force"] + faith.force.youth$coefficients["force:exposureself.interned"], faith.force.youth$coefficients["force"],
                                      leader.force.youth$coefficients["force"] + leader.force.youth$coefficients["force:exposureself.interned"], leader.force.youth$coefficients["force"], 
                                      polad.force.youth$coefficients["force"] + polad.force.youth$coefficients["force:exposureself.interned"], polad.force.youth$coefficients["force"]),
                          ses = c(sqrt(interest.force.youth.ses[[2]]["force", "force"] + interest.force.youth.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * interest.force.youth.ses[[2]]["force", "force:exposureself.interned"]) , 
                                  sqrt(interest.force.youth.ses[[2]]["force", "force"]),
                                  sqrt(faith.force.youth.ses[[2]]["force", "force"] + faith.force.youth.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * faith.force.youth.ses[[2]]["force", "force:exposureself.interned"]) , 
                                  sqrt(faith.force.youth.ses[[2]]["force", "force"]), 
                                  sqrt(leader.force.youth.ses[[2]]["force", "force"] + leader.force.youth.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * leader.force.youth.ses[[2]]["force", "force:exposureself.interned"]) , 
                                  sqrt(leader.force.youth.ses[[2]]["force", "force"]),
                                  sqrt(polad.force.youth.ses[[2]]["force", "force"] + polad.force.youth.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * polad.force.youth.ses[[2]]["force", "force:exposureself.interned"]) , 
                                  sqrt(polad.force.youth.ses[[2]]["force", "force"])))
  force.youth.res$ub <- force.youth.res$results + qnorm(.975) * force.youth.res$ses
  force.youth.res$lb <- force.youth.res$results - qnorm(.975) * force.youth.res$ses
  force.youth.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  force.youth.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  force.youth.res$labels <- factor(force.youth.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  force.youth.res$groups <- factor(force.youth.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  force.la.res <- data.frame(results = c(interest.force.la$coefficients["force"] + interest.force.la$coefficients["force:exposureself.interned"], interest.force.la$coefficients["force"],
                                            faith.force.la$coefficients["force"] + faith.force.la$coefficients["force:exposureself.interned"], faith.force.la$coefficients["force"],
                                            leader.force.la$coefficients["force"] + leader.force.la$coefficients["force:exposureself.interned"], leader.force.la$coefficients["force"], 
                                            polad.force.la$coefficients["force"] + polad.force.la$coefficients["force:exposureself.interned"], polad.force.la$coefficients["force"]),
                                ses = c(sqrt(interest.force.la.ses[[2]]["force", "force"] + interest.force.la.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * interest.force.la.ses[[2]]["force", "force:exposureself.interned"]) , 
                                        sqrt(interest.force.la.ses[[2]]["force", "force"]),
                                        sqrt(faith.force.la.ses[[2]]["force", "force"] + faith.force.la.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * faith.force.la.ses[[2]]["force", "force:exposureself.interned"]) , 
                                        sqrt(faith.force.la.ses[[2]]["force", "force"]), 
                                        sqrt(leader.force.la.ses[[2]]["force", "force"] + leader.force.la.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * leader.force.la.ses[[2]]["force", "force:exposureself.interned"]) , 
                                        sqrt(leader.force.la.ses[[2]]["force", "force"]),
                                        sqrt(polad.force.la.ses[[2]]["force", "force"] + polad.force.la.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * polad.force.la.ses[[2]]["force", "force:exposureself.interned"]) , 
                                        sqrt(polad.force.la.ses[[2]]["force", "force"])))
  force.la.res$ub <- force.la.res$results + qnorm(.975) * force.la.res$ses
  force.la.res$lb <- force.la.res$results - qnorm(.975) * force.la.res$ses
  force.la.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  force.la.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  force.la.res$labels <- factor(force.la.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  force.la.res$groups <- factor(force.la.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  
  force.birthplace.res <- data.frame(results = c(interest.force.birthplace$coefficients["force"] + interest.force.birthplace$coefficients["force:exposureself.interned"], interest.force.birthplace$coefficients["force"],
                                      faith.force.birthplace$coefficients["force"] + faith.force.birthplace$coefficients["force:exposureself.interned"], faith.force.birthplace$coefficients["force"],
                                      leader.force.birthplace$coefficients["force"] + leader.force.birthplace$coefficients["force:exposureself.interned"], leader.force.birthplace$coefficients["force"], 
                                      polad.force.birthplace$coefficients["force"] + polad.force.birthplace$coefficients["force:exposureself.interned"], polad.force.birthplace$coefficients["force"]),
                          ses = c(sqrt(interest.force.birthplace.ses[[2]]["force", "force"] + interest.force.birthplace.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * interest.force.birthplace.ses[[2]]["force", "force:exposureself.interned"]) , 
                                  sqrt(interest.force.birthplace.ses[[2]]["force", "force"]),
                                  sqrt(faith.force.birthplace.ses[[2]]["force", "force"] + faith.force.birthplace.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * faith.force.birthplace.ses[[2]]["force", "force:exposureself.interned"]) , 
                                  sqrt(faith.force.birthplace.ses[[2]]["force", "force"]), 
                                  sqrt(leader.force.birthplace.ses[[2]]["force", "force"] + leader.force.birthplace.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * leader.force.birthplace.ses[[2]]["force", "force:exposureself.interned"]) , 
                                  sqrt(leader.force.birthplace.ses[[2]]["force", "force"]),
                                  sqrt(polad.force.birthplace.ses[[2]]["force", "force"] + polad.force.birthplace.ses[[2]]["force:exposureself.interned", "force:exposureself.interned"] + 2 * polad.force.birthplace.ses[[2]]["force", "force:exposureself.interned"]) , 
                                  sqrt(polad.force.birthplace.ses[[2]]["force", "force"])))
  force.birthplace.res$ub <- force.birthplace.res$results + qnorm(.975) * force.birthplace.res$ses
  force.birthplace.res$lb <- force.birthplace.res$results - qnorm(.975) * force.birthplace.res$ses
  force.birthplace.res$labels <- c(rep("Political \n Interest", 2), rep("Political \n Distrust", 2), rep("Leadership \n Approach", 2), rep("Political \n Advice", 2))
  force.birthplace.res$groups <- rep(c("Direct Exposure", "Family-Only Exposure"), 4)
  force.birthplace.res$labels <- factor(force.birthplace.res$labels, levels = c("Leadership \n Approach", "Political \n Advice", "Political \n Distrust", "Political \n Interest"), ordered = T)
  force.birthplace.res$groups <- factor(force.birthplace.res$groups, levels = c("Family-Only Exposure", "Direct Exposure"), ordered = T)
  
  
  pdf("ms_fig9b.pdf")
    ggplot(data = force.res, aes(x = labels, y = results, color = groups)) + 
      geom_hline(yintercept = 0, alpha = .5) + 
      geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
      geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
      theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
      labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Use of Military Force While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig16c.pdf")
    ggplot(data = force.youth.res, aes(x = labels, y = results, color = groups)) + 
      geom_hline(yintercept = 0, alpha = .5) + 
      geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
      geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
      theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
      labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Use of Military Force While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig15c.pdf")
    ggplot(data = force.la.res, aes(x = labels, y = results, color = groups)) + 
      geom_hline(yintercept = 0, alpha = .5) + 
      geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
      geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
      theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
      labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Use of Military Force While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  pdf("si_fig14c.pdf")
  ggplot(data = force.birthplace.res, aes(x = labels, y = results, color = groups)) + 
    geom_hline(yintercept = 0, alpha = .5) + 
    geom_pointrange(aes(ymin = lb, ymax = ub), size = .5, position = position_dodge(.5)) + 
    geom_errorbar(aes(ymin = results - 1.405 * ses, ymax = results + 1.405 * ses), size = 1, width = 0, position = position_dodge(.5)) + 
    theme_minimal() + scale_color_grey(start = .6, end = .1) + theme(plot.title = element_text(hjust = 0.5)) +
    labs(x = "Outcome", y = "Coefficient Estimate", title = "Effect of Witnessing Use of Military Force While Interned", color = "Group") + coord_flip() +  guides(color = guide_legend(reverse = TRUE))
  dev.off()
  
  
# Age Interaction Plots
  
# Interest  
  interest.age.plots.violence <- interplot(m = interest.violence.ageinter, var1 = "violence", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Interest in Politics", title = "Violence x Age") 

  interest.age.plots.strikes <- interplot(m = interest.strikes.ageinter, var1 = "strikes", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Interest in Politics", title = "Demonstrations x Age") 
  
  interest.age.plots.force <- interplot(m = interest.force.ageinter, var1 = "force", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Interest in Politics", title = "Force x Age") 
  
  interest.age.plots.militarism <- interplot(m = interest.militarism.ageinter, var1 = "militarism", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Interest in Politics", title = "Militarism x Age") 
  
  pdf("si_fig19a.pdf")
    grid.arrange(interest.age.plots.violence, interest.age.plots.strikes, interest.age.plots.force, interest.age.plots.militarism)
  dev.off()
  
  
# Political Distrust
  faith.age.plots.violence <- interplot(m = faith.violence.ageinter, var1 = "violence", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Political Distrust", title = "Violence x Age") 
  
  faith.age.plots.strikes <- interplot(m = faith.strikes.ageinter, var1 = "strikes", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Political Distrust", title = "Demonstrations x Age") 
  
  faith.age.plots.force <- interplot(m = faith.force.ageinter, var1 = "force", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Political Distrust", title = "Force x Age") 
  
  faith.age.plots.militarism <- interplot(m = faith.militarism.ageinter, var1 = "militarism", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Political Distrust", title = "Militarism x Age")   
  
  pdf("si_fig19b.pdf")
    grid.arrange(faith.age.plots.violence, faith.age.plots.strikes, faith.age.plots.force, faith.age.plots.militarism)
  dev.off()
  
  
# Leadership
  leader.age.plots.violence <- interplot(m = leader.violence.ageinter, var1 = "violence", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Leadership Approach", title = "Violence x Age") 
  
  leader.age.plots.strikes <- interplot(m = leader.strikes.ageinter, var1 = "strikes", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Leadership Approach", title = "Demonstrations x Age") 
  
  leader.age.plots.force <- interplot(m = leader.force.ageinter, var1 = "force", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Leadership Approach", title = "Force x Age") 
  
  leader.age.plots.militarism <- interplot(m = leader.militarism.ageinter, var1 = "militarism", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Leadership Approach", title = "Militarism x Age")   
  
  pdf("si_fig19c.pdf")
    grid.arrange(leader.age.plots.violence, leader.age.plots.strikes, leader.age.plots.force, leader.age.plots.militarism)
  dev.off()
  
# Political Advice
  polad.age.plots.violence <- interplot(m = polad.violence.ageinter, var1 = "violence", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Political Advice", title = "Violence x Age") 
  
  polad.age.plots.strikes <- interplot(m = polad.strikes.ageinter, var1 = "strikes", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Political Advice", title = "Demonstrations x Age") 
  
  polad.age.plots.force <- interplot(m = polad.force.ageinter, var1 = "force", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Political Advice", title = "Force x Age") 
  
  polad.age.plots.militarism <- interplot(m = polad.militarism.ageinter, var1 = "militarism", var2 = "age") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Age", y = "Political Advice", title = "Militarism x Age")   
  
  pdf("si_fig19d.pdf")
    grid.arrange(polad.age.plots.violence, polad.age.plots.strikes, polad.age.plots.force, polad.age.plots.militarism)
  dev.off()  
  
  # Interest  
  interest.minors.plots.violence <- interplot(m = interest.violence.minors, var1 = "violence", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Interest in Politics", title = "Violence x Minor") 
  
  interest.minors.plots.strikes <- interplot(m = interest.strikes.minors, var1 = "strikes", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Interest in Politics", title = "Demonstrations x Minor") 
  
  interest.minors.plots.force <- interplot(m = interest.force.minors, var1 = "force", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Interest in Politics", title = "Force x Minor") 
  
  interest.minors.plots.militarism <- interplot(m = interest.militarism.minors, var1 = "militarism", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Interest in Politics", title = "Militarism x Minor") 
  
  pdf("si_fig20a.pdf")
    grid.arrange(interest.minors.plots.violence, interest.minors.plots.strikes, interest.minors.plots.force, interest.minors.plots.militarism)
  dev.off()
  
  # Political Distrust
  faith.minors.plots.violence <- interplot(m = faith.violence.minors, var1 = "violence", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Political Distrust", title = "Violence x Minor") 
  
  faith.minors.plots.strikes <- interplot(m = faith.strikes.minors, var1 = "strikes", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Political Distrust", title = "Demonstrations x Minor") 
  
  faith.minors.plots.force <- interplot(m = faith.force.minors, var1 = "force", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Political Distrust", title = "Force x Minor") 
  
  faith.minors.plots.militarism <- interplot(m = faith.militarism.minors, var1 = "militarism", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Political Distrust", title = "Militarism x Minor")   
  
  pdf("si_fig20b.pdf")
    grid.arrange(faith.minors.plots.violence, faith.minors.plots.strikes, faith.minors.plots.force, faith.minors.plots.militarism)
  dev.off()
  
  
  # Leadership
  leader.minors.plots.violence <- interplot(m = leader.violence.minors, var1 = "violence", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Leadership Approach", title = "Violence x Minor") 
  
  leader.minors.plots.strikes <- interplot(m = leader.strikes.minors, var1 = "strikes", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Leadership Approach", title = "Demonstrations x Minor") 
  
  leader.minors.plots.force <- interplot(m = leader.force.minors, var1 = "force", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Leadership Approach", title = "Force x Minor") 
  
  leader.minors.plots.militarism <- interplot(m = leader.militarism.minors, var1 = "militarism", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Leadership Approach", title = "Militarism x Minor")   
  
  pdf("si_fig20c.pdf")
    grid.arrange(leader.minors.plots.violence, leader.minors.plots.strikes, leader.minors.plots.force, leader.minors.plots.militarism)
  dev.off()
  
  # Political Advice
  polad.minors.plots.violence <- interplot(m = polad.violence.minors, var1 = "violence", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Political Advice", title = "Violence x Minor") 
  
  polad.minors.plots.strikes <- interplot(m = polad.strikes.minors, var1 = "strikes", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Political Advice", title = "Demonstrations x Minor") 
  
  polad.minors.plots.force <- interplot(m = polad.force.minors, var1 = "force", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Political Advice", title = "Force x Minor") 
  
  polad.minors.plots.militarism <- interplot(m = polad.militarism.minors, var1 = "militarism", var2 = "minordummy") + theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5)) +
    geom_hline(yintercept = 0, alpha = .5, linetype = "dashed") + 
    labs(x = "Minor", y = "Political Advice", title = "Militarism x Minor")   
  
  pdf("si_fig20d.pdf")
    grid.arrange(polad.minors.plots.violence, polad.minors.plots.strikes, polad.minors.plots.force, polad.minors.plots.militarism)
  dev.off()  
  
  
# Export Results as Appendix Table
  stargazer(interest.violence, faith.violence, polad.violence, leader.violence,
            interest.strikes, faith.strikes, polad.strikes, leader.strikes,
            interest.force, faith.force, polad.force, leader.force,
            interest.militarism, faith.militarism, polad.militarism, leader.militarism,
            se = list(sqrt(diag(interest.violence.ses[[2]])), sqrt(diag(faith.violence.ses[[2]])), sqrt(diag(polad.violence.ses[[2]])), sqrt(diag(leader.violence.ses[[2]])),
                      sqrt(diag(interest.strikes.ses[[2]])), sqrt(diag(faith.strikes.ses[[2]])), sqrt(diag(polad.strikes.ses[[2]])), sqrt(diag(leader.strikes.ses[[2]])),
                      sqrt(diag(interest.force.ses[[2]])), sqrt(diag(faith.force.ses[[2]])), sqrt(diag(polad.force.ses[[2]])), sqrt(diag(leader.force.ses[[2]])),
                      sqrt(diag(interest.militarism.ses[[2]])), sqrt(diag(faith.militarism.ses[[2]])), sqrt(diag(polad.militarism.ses[[2]])), sqrt(diag(leader.militarism.ses[[2]]))),
            dep.var.labels = rep(c("Interest", "Distrust", "Advice", "Leadership"), 4),
            font.size = "tiny",
            no.space = T,
            single.row = T,
            digits = 2,
            column.sep.width = "-15pt",
            title = "Effects of Internment Experience", 
            initial.zero = FALSE,
            keep.stat = c("n", "rsq"),
            label = "tab:fullrestab", table.placement = "H",
            out = "si_table11.tex",
            add.lines = list(c("Generations", "I, N, S", "N, S", "N, S", "N, S", "I, N, S", "N, S", "N, S", "N, S", "I, N, S", "N, S", "N, S", "N, S", "I, N, S", "N, S", "N, S", "N, S"),
                             c("Generation Fixed Effects", 
                               "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark"),
                             c("Pre-Internment Location Fixed Effects", 
                               "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark")),
            covariate.labels = c("Violence", "Strikes", "Force", "Militarism","Interned","Gender (Male)", "Age", "Violence x Interned", "Demonstrations x Interned", "Force x Interned", "Militarism x Interned"),
            keep = c("violence", "strikes", "force", "militarism", "violence:exposureself.interned", "strikes:exposureself.interned", "force:exposureself.interned", "militarism:exposureself.interned", "exposureself.interned", "age", "genderMale"),
            notes = c("Note: Standard errors clustered by internment camp to which respondents are assigned."))     

  
################################################################################################################################################################################################################################################
# Robustness Checks: Covariate Balance
################################################################################################################################################################################################################################################
# Add Family Size (prior to internment)
# Get number of family members by ID who are 21 + (would have been alive in 41 given that survey is 1962)
famsize <- data.frame(table(full$respid[full$age > 21]))
colnames(famsize) <- c("respid", "pt.fam.members")
full <- merge(full, famsize, by = "respid", all.x = T)

# Balance Test
# Simplest form: Let's Regress Camp in Multinomial Logit form on Gender, Age, Family Size, Pre-Camp Location
    camp.assignment <- multinom(camp ~ age + gender + as.numeric(precamp), data = full[!(full$camp %in% c("Not evacuated or interned", "No Data", "Livingston, LA", "Fort Sill, OK", "Missoula, MT", "Crystal City, TX", "Bismark, ND", "Lordsburg, NM", "Santa Fe, NM")) & full$treatmentstatus %in% c("self.only", "self.fam"),])  
    z <- data.frame(summary(camp.assignment)$coefficients[ , -1] / summary(camp.assignment)$standard.errors[ , -1])
    colnames(z) <-  c("Age", "Gender (Male)", "Pre-Camp Location") 
    rownames(z) <- c("Jerome", "Heart Mountain", "Minidoka", "Manzanar", "Rohwer", "Tule Lake", "Poston", "Gila River", "Topaz")
     
    stargazer(z, summary = F, font.size = "tiny", rownames = T, digits = 2,
              title = "Covariate Balance Across Internment Camps", no.space = T,
              dep.var.caption = "",
              dep.var.labels = "",
              initial.zero = FALSE,
              notes = c("Note: Reference category is Granada (Amache)"),
              out = "ms_table4.tex",
              label = "tab:exomult", table.placement = "H")   


# KS Tests for Age, Gender, Employment
  # Get all possible unique pairs of camps
  
  # Unique camps
    camp.pairs <- unique(full$camp[!(full$camp %in% c("Not evacuated or interned", "No Data", "Livingston, LA", "Fort Sill, OK", "Missoula, MT", "Crystal City, TX", "Bismark, ND", "Lordsburg, NM", "Santa Fe, NM")) & !(is.na(full$camp)) & full$treatmentstatus %in% c("self.only", "self.fam")])
    
  # All pairs  
    camp.pairs <- expand.grid(camp.pairs, camp.pairs)
    
  # Drop same camp
    camp.pairs <- camp.pairs[camp.pairs$Var1 != camp.pairs$Var2,] 
    
  # Drop duplicates
    ks.camps <- camp.pairs[1,]
    
    for(i in 2:nrow(camp.pairs)){
      if(!(camp.pairs$Var1[i] %in% camp.pairs$Var2[1:i])){
        ks.camps <- rbind(ks.camps, camp.pairs[i,])
      }
    }
  # Roll through age to get a series of p-values on K.S. tests by camp pair
    ks.age <- NULL
    for(i in 1:nrow(ks.camps)){
      ks.age[i] <- ks.test(x = full$age[full$camp == ks.camps$Var1[i] & full$treatmentstatus %in% c("self.only", "self.fam")], y = full$age[full$camp == ks.camps$Var2[i] & full$treatmentstatus %in% c("self.only", "self.fam")])$p.value
    }
      
  # Same thing for jobs, except that I'm going to separate off Nisei since that's the only place jobs is personal job, not head of household
    ks.ptjobs.n <- NULL
    for(i in 1:nrow(ks.camps)){
      ks.ptjobs.n[i] <- ks.test(x = full$n.ptjob[full$camp == ks.camps$Var1[i]& full$treatmentstatus %in% c("self.only", "self.fam")], y = full$n.ptjob[full$camp == ks.camps$Var2[i]& full$treatmentstatus %in% c("self.only", "self.fam")])$p.value
    }
    
    age.ks <- data.frame(cbind(ks.camps, ks.age))
    colnames(age.ks) <- c("Camp 1", "Camp 2", "p-value")
    
    stargazer(age.ks, summary = F, font.size = "tiny", rownames = F, digits = 2,
              title = "Kolmogorov-Smirnov Test Results for Age", no.space = T,
              dep.var.caption = "", 
              dep.var.labels = "",
              out = "si_table12.tex",
              label = "tab:ksage", table.placement = "H")
    
    ptjob.ks <- data.frame(cbind(ks.camps,  ks.ptjobs.n))
    colnames(ptjob.ks) <- c("Camp 1", "Camp 2", "p-value")
    
    stargazer(ptjob.ks, summary = F, font.size = "tiny", rownames = F, digits = 2,
              title = "Kolmogorov-Smirnov Test for Employment Sector", no.space = T,
              dep.var.caption = "", 
              dep.var.labels = "",
              out = "si_table13.tex",
              label = "tab:ksptjob", table.placement = "H")

# Standardized differences in means by gender
  genderdiff <- data.frame(diffmean = rep(NA, nrow(ks.camps)), se =  rep(NA, nrow(ks.camps)), t = rep(NA, nrow(ks.camps)), p = rep(NA, nrow(ks.camps)), lb = rep(NA, nrow(ks.camps)), ub = rep(NA, nrow(ks.camps)))
  
  for(i in 1:nrow(ks.camps)){
    x <- full[full$camp == ks.camps$Var1[i] & full$treatmentstatus %in% c("self.only", "self.fam"),]
    y <- full[full$camp == ks.camps$Var2[i] & full$treatmentstatus %in% c("self.only", "self.fam"),]
    p <- mean(full$gender[full$camp %in% c(as.character(ks.camps$Var1[i]), as.character(ks.camps$Var2[i]))] == 0, na.rm = T)
    genderdiff$diffmean[i] <- mean(y$gender == 0, na.rm = T) - mean(x$gender == 0, na.rm = T)
    genderdiff$se[i] <-  sqrt(p * (1 - p) / (1 / length(!is.na(y$gender)) + 1 / length(!is.na(x$gender))))
    genderdiff$t[i] <- genderdiff$diffmean[i] / genderdiff$se[i]
    genderdiff$lb[i] <- genderdiff$t[i] - qt(0.975, df = min(length(!is.na(x$gender)) - 1, length(!is.na(y$gender)) - 1)) * genderdiff$se[i]
    genderdiff$ub[i] <- genderdiff$t[i] + qt(0.975, df = min(length(!is.na(x$gender)) - 1, length(!is.na(y$gender)) - 1)) * genderdiff$se[i]
    genderdiff$p[i] <- ifelse(genderdiff$t[i] > 0, 2 * (1 - pt(genderdiff$t[i], df =  min(length(!is.na(x$gender)) - 1, length(!is.na(y$gender)) - 1))),
                               2 * (pt(genderdiff$t[i], df =  min(length(!is.na(x$gender)) - 1, length(!is.na(y$gender)) - 1))))
  }
  
  genderdiff <- data.frame(cbind(ks.camps, genderdiff))
  
  colnames(genderdiff) <- c("Camp 1", "Camp 2","Difference in Means", "Standard Error", "T", "p-value", "Lower Bound", "Upper Bound")
  
  stargazer(genderdiff, summary = F, font.size = "tiny", rownames = F, digits = 2,
            title = "T-Tests for Differences in Gender Proportion by Camp", no.space = T,
            dep.var.caption = "", 
            dep.var.labels = "",
            out = "si_table14.tex",
            label = "tab:gendertest", table.placement = "H")

  
################################################################################################################################################################################################################################################
# Robustness Checks: Alternative Mechanisms
################################################################################################################################################################################################################################################
# Load 1940 Census Data
  c40 <- read_dta("1940b.dta")
  c40 <- c40[, c("fips", "name", "totpop", "mtot", "ftot", "wmtot", "wftot")]
  c40$pcw40 <- (c40$wmtot + c40$wftot) / c40$totpop * 100

# Merge To Counties with Camp Locations
  data(state.fips)
  c40$statefip <- ifelse(nchar(c40$fips) == 4, substr(c40$fips, 1, 1), ifelse(nchar(c40$fips) == 5, substr(c40$fips, 1, 2), NA))
  c40 <- merge(c40, state.fips, by.x = "statefip", by.y = "fips")
  c40$polyname <- toupper(c40$polyname)
  camps <- merge(camps, c40, by.x = c("county", "state"), by.y = c("name", "polyname"))

# Load Historical Political Data from ICPSR
  votes <- read.dta("08611-0001-Data.dta")
  colnames(votes) <- attr(votes, "var.labels")

  votes$state <- trimws(splitit(votes$`ICPSR STATE CODE`, "\\(", 1), which = "right")
  colnames(votes)[2] <- "county"

  votes2 <- merge(camps, votes, by = c("state", "county"), all.x = T)
  votes2 <- votes2[ , colnames(votes2)[grepl("state", colnames(votes2)) == 1 | grepl("county", colnames(votes2)) == 1 | 
                                       grepl("camp", colnames(votes2)) == 1 | grepl("1940", colnames(votes2)) == 1]]
  votes2$presdem40 <- votes2$`1940 PRES TTL VOTE` * votes2$`1940 PRES DEM          PERCENT` / 
    (votes2$`1940 PRES TTL VOTE` * votes2$`1940 PRES DEM          PERCENT` + votes2$`1940 PRES TTL VOTE` * votes2$`1940 PRES REP          PERCENT`)

  votes2 <- votes2[ , c("camp", "presdem40")]

  full <- merge(full, votes2, by = "camp", all.x = T)
  full <- merge(full, camps[,c("camp", "pcw40")], by = "camp", all.x = T)
  

# Robustness to Including % White and Local Vote Share
    # Violence  
      # Complete Data for Standard Errors - include % white + FDR vote share
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("interestcodes", "pcw40", "presdem40", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        interest.violence.robust <- lm(interestcodes ~ violence * exposure + pcw40 + presdem40 + precamp + gender + age + generation, data = comp)  
      
      # Clustered Standard Errors 
        interest.violence.robust.ses  <- cl(comp, interest.violence.robust, comp$camp)
    
    # Force
      # Complete Data for Standard Errors - include % white + FDR vote share
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("interestcodes", "pcw40", "presdem40", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        interest.force.robust <- lm(interestcodes ~ force * exposure + pcw40 + presdem40 + precamp + gender + age + generation, data = comp)  
        
      # Clustered Standard Errors 
        interest.force.robust.ses  <- cl(comp, interest.force.robust, comp$camp)
        
    # Strikes
      # Complete Data for Standard Errors - include % white + FDR vote share
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("interestcodes", "pcw40", "presdem40", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        interest.strikes.robust <- lm(interestcodes ~ strikes * exposure + pcw40 + presdem40 + precamp + gender + age + generation, data = comp)  
        
      # Clustered Standard Errors 
        interest.strikes.robust.ses  <- cl(comp, interest.strikes.robust, comp$camp)
        
    # Militarism
      # Complete Data for Standard Errors - include % white + FDR vote share
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("interestcodes", "pcw40", "presdem40", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        interest.militarism.robust <- lm(interestcodes ~ militarism * exposure + pcw40 + presdem40 + precamp + gender + age + generation, data = comp)  
        
      # Clustered Standard Errors 
        interest.militarism.robust.ses  <- cl(comp, interest.militarism.robust, comp$camp)
        
  # Faith in Government - just strikes significant here (+)
    # Violence  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("govhelp", "pcw40", "presdem40", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        faith.violence.robust <- lm(govhelp ~ violence * exposure + pcw40 + presdem40 + precamp + age + gender + generation, data = comp)
        
      # Clustered Standard Errors 
        faith.violence.robust.ses  <- cl(comp, faith.violence.robust, comp$camp)
        
    # Force  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("govhelp", "pcw40", "presdem40", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        faith.force.robust <- lm(govhelp ~ force * exposure + pcw40 + presdem40 + precamp + age + gender + generation, data = comp)
        
      # Clustered Standard Errors 
        faith.force.robust.ses  <- cl(comp, faith.force.robust, comp$camp)
        
    # Strikes  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("govhelp", "pcw40", "presdem40", "strikes", "precamp", "age", "gender", "camp", "generation","exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        faith.strikes.robust <- lm(govhelp ~ strikes * exposure + pcw40 + presdem40 + precamp + age + gender + generation, data = comp)
        
      # Clustered Standard Errors 
        faith.strikes.robust.ses  <- cl(comp, faith.strikes.robust, comp$camp)
      
    # Militarism  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("govhelp", "pcw40", "presdem40", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        faith.militarism.robust <- lm(govhelp ~ militarism * exposure + pcw40 + presdem40 + precamp + age + gender + generation, data = comp)
        
      # Clustered Standard Errors 
        faith.militarism.robust.ses  <- cl(comp, faith.militarism.robust, comp$camp)
        
  # Leadership - just violence and strikes here (both +)
    # Violence  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("leader", "pcw40", "presdem40", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        leader.violence.robust <- lm(leader ~ violence * exposure + pcw40 + presdem40 + precamp + gender + age + generation, data = comp)  
        
      # Clustered Standard Errors 
        leader.violence.robust.ses  <- cl(comp, leader.violence.robust, comp$camp)
        
    # Force  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("leader", "pcw40", "presdem40", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        leader.force.robust <- lm(leader ~ force * exposure + pcw40 + presdem40 + precamp + gender + age + generation, data = comp)  
        
      # Clustered Standard Errors 
        leader.force.robust.ses  <- cl(comp, leader.force.robust, comp$camp)
        
    # Strikes  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("leader", "pcw40", "presdem40", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        leader.strikes.robust <- lm(leader ~ strikes * exposure + pcw40 + presdem40 + precamp + gender + age + generation, data = comp)  
        
      # Clustered Standard Errors 
        leader.strikes.robust.ses  <- cl(comp, leader.strikes.robust, comp$camp)
        
    # Militarism  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("leader", "pcw40", "presdem40", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        leader.militarism.robust <- lm(leader ~ militarism * exposure + pcw40 + presdem40 + precamp + gender + age + generation, data = comp)  
        
      # Clustered Standard Errors 
        leader.militarism.robust.ses  <- cl(comp, leader.militarism.robust, comp$camp)
        
  # Political Advice - just militarism (+) and strikes (-) here         
    # Violence  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("poladvc", "pcw40", "presdem40", "violence", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        polad.violence.robust <- lm(poladvc ~ violence * exposure + pcw40 + presdem40 + precamp + age + gender + generation, data = comp)
        
      # Clustered Standard Errors 
        polad.violence.robust.ses  <- cl(comp, polad.violence.robust, comp$camp)
        
    # Force  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("poladvc", "pcw40", "presdem40", "force", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        polad.force.robust <- lm(poladvc ~ force * exposure + pcw40 + presdem40 + precamp + age + gender + generation, data = comp)
        
      # Clustered Standard Errors 
        polad.force.robust.ses  <- cl(comp, polad.force.robust, comp$camp)
        
    # Strikes  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("poladvc", "pcw40", "presdem40", "strikes", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
      
      # Model  
        polad.strikes.robust <- lm(poladvc ~ strikes * exposure + pcw40 + presdem40 + precamp + age + gender + generation, data = comp)
      
      # Clustered Standard Errors 
        polad.strikes.robust.ses  <- cl(comp, polad.strikes.robust, comp$camp)
      
    # Militarism  
      # Complete Data for Standard Errors
        comp <- full[full$exposure %in% c("fam.only", "self.interned"), c("poladvc", "pcw40", "presdem40", "militarism", "precamp", "age", "gender", "camp", "generation", "exposure")]
        comp <- comp[complete.cases(comp),]
        
      # Model  
        polad.militarism.robust <- lm(poladvc ~ militarism * exposure + pcw40 + presdem40 + precamp + age + gender + generation, data = comp)
        
      # Clustered Standard Errors 
        polad.militarism.robust.ses  <- cl(comp, polad.militarism.robust, comp$camp)
        
        stargazer(interest.violence.robust, faith.violence.robust, polad.violence.robust, leader.violence.robust,
                  interest.strikes.robust, faith.strikes.robust, polad.strikes.robust, leader.strikes.robust,
                  interest.force.robust, faith.force.robust, polad.force.robust, leader.force.robust,
                  interest.militarism.robust, faith.militarism.robust, polad.militarism.robust, leader.militarism.robust,
                  se = list(sqrt(diag(interest.violence.robust.ses[[2]])), sqrt(diag(faith.violence.robust.ses[[2]])), sqrt(diag(polad.violence.robust.ses[[2]])), sqrt(diag(leader.violence.robust.ses[[2]])),
                            sqrt(diag(interest.strikes.robust.ses[[2]])), sqrt(diag(faith.strikes.robust.ses[[2]])), sqrt(diag(polad.strikes.robust.ses[[2]])), sqrt(diag(leader.strikes.robust.ses[[2]])),
                            sqrt(diag(interest.force.robust.ses[[2]])), sqrt(diag(faith.force.robust.ses[[2]])), sqrt(diag(polad.force.robust.ses[[2]])), sqrt(diag(leader.force.robust.ses[[2]])),
                            sqrt(diag(interest.militarism.robust.ses[[2]])), sqrt(diag(faith.militarism.robust.ses[[2]])), sqrt(diag(polad.militarism.robust.ses[[2]])), sqrt(diag(leader.militarism.robust.ses[[2]]))),
                  dep.var.labels = rep(c("Interest", "Distrust", "Advice", "Leadership"), 4),
                  font.size = "tiny",
                  no.space = T,
                  single.row = T,
                  digits = 2,
                  column.sep.width = "-13.5pt",
                  title = "Effects of Internment Experience, Controlling for Surrounding Camp Environment", 
                  initial.zero = FALSE,
                  keep.stat = c("n", "rsq"),
                  label = "tab:fullrestabrobust", table.placement = "H",
                  out = "si_table16.tex",
                  add.lines = list(c("Generations", "I, N, S", "N, S", "N, S", "N, S", "I, N, S", "N, S", "N, S", "N, S", "I, N, S", "N, S", "N, S", "N, S", "I, N, S", "N, S", "N, S", "N, S"),
                                   c("Generation FE", 
                                     "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark"),
                                   c("Location FE", 
                                     "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark", "\\checkmark")),
                  covariate.labels = c("Violence", "Strikes", "Force", "Militarism","Interned", "Pct. White 1940", "FDR 1940 Vote Share","Gender (Male)", "Age", "Violence x Interned", "Demonstrations x Interned", "Force x Interned", "Militarism x Interned"),
                  keep = c("violence", "strikes", "force", "militarism", "violence:exposureself.interned", "strikes:exposureself.interned", "force:exposureself.interned", "militarism:exposureself.interned", "exposureself.interned", "age", "genderMale", "pcw40", "presdem40"),
                  notes = c("Note: Standard errors clustered by internment camp to which respondents are assigned. Location FE refer to pre-internment location fixed effects."))       